function [dcdt, algvars] = MnPfunc(t,c,params)
% Assign names for parameter values and state variables
[p,v,b0,j,k1f,k1r,k2] = params{:};
MnP = c(1);
Bac = c(2);
Nyl = c(3);
NylMnP= c(4);

% Conservation equations
FreeMnP = MnP - NylMnP;                                                    % [unit] free MnP concentration
FreeNyl = Nyl - NylMnP;                                                    % [unit] free MnP concentration

% Reaction rates
react1 = p*b0/v- FreeMnP*j;                                                % production and elimination of MnP
react2 = Bac;                                                              % Assume bacteria reaches saturation/confluence
react3 = k1f*FreeNyl*FreeMnP-k1r*FreeNyl*FreeMnP;                          % [unit/s] Nylon + MnP <-> MnPNyl
react4 = k2*NylMnP;                                                        % [unit/s] NylMnP <-> product + MnP
react5 = Nyl*j/v;                                                          % Nylon constantly flow into the chamber

% Differential equations;
dMnP = react1;                                                             % [unit/min] MnP rate of change
dBac = 0;                                                                  % [unit/min] Bacteria rate of change 
dNyl = react5-react3;                                                      % [unit/min] Nylon rate of change
dNylMnP = react3 - react4;                                                 % [unit/min] Nylon-MnP rate of change
dcdt = [dMnP;dBac;dNyl;dNylMnP];                                           % Reassemble differential equations
algvars = [FreeMnP,react1,react2,react3,react4,react5];