MATLAB: Implementing a code from Berkley Madonna into Matlab

berkely madonnadifferential equationsmitotic oscillator

I am trying to implement a code from Berkley Madonna into Matlab. I want to carry out a simulation and produce the results graphically by using ode solvers directly. Here is this Berkley Madonna code I am trying to incorporate:
Listing of program in Berkeley Madonna notation:
DT=1e-4
STOPTIME=100
d/dt (C) = Synthesis – Degradation
INIT C = 0.01
Synthesis = 0.025
Degradation = vd*X*(C/(Kd+C)) – kdd*C
d/dt (M) = Phosphatase1 – Kinase1
INIT M = 0.01
Phosphatase1 = VM1*(C/(Kc+C))*((1-M)/(K1+(1-M)))
Kinase1 = V2*(M/(K2+M))
d/dt (X) = Phosphatase2 – Kinase2
INIT X = 0.01
Phosphatase2 = M*VM3*((1-X)/(K3+(1-X)))
Kinase2 = V4*(X/(K4+X))
K1 = 0.005
K2 = 0.005
K3 = 0.005
K4 = 0.005
Kc = 0.5
Kd = 0.02
kdd = 0.01
V2 = 1.5
V4 = 0.5
vd = 0.25
VM1 = 3
VM3 = 1.
I have first two questions about this: First is the equation set up given by Berkley Madonna in correct use for matlab. Second question….I am unable to use the ode solvers correctly. It seems I somehow have re written over the ode solvers that are built into Matlab. Is there a way to correct this without uninstalling MATLAB and re installing it.

Best Answer

Just for fun (and because I’m interested in biochemistry), I decided to code this as best I could. It gives interesting results, but certainly not out of character for some biochemical systems. (I copied the equations directly from your code, pasted them into the function, and did a ‘Search - Replace’ with the Editor.)
K1 = 0.005;
K2 = 0.005;
K3 = 0.005;
K4 = 0.005;
Kc = 0.5;
Kd = 0.02;
kdd = 0.01;
V2 = 1.5;
V4 = 0.5;
vd = 0.25;
VM1 = 3;
VM3 = 1.;
Synth = 0.025;
% Original Equations: 
% PhsKnsKnt = [Synth-vd*X*(C/(Kd+C)) - kdd*C;    M*VM3*((1-X)/(K3+(1-X))) - V4*(X/(K4+X));    VM1*(C/(Kc+C))*((1-M)/(K1+(1-M))) - V2*(M/(K2+M))];
% Anonymous Function with substitutions: 
PhsKnsKnt = @(T, P)  [Synth - vd*P(2)*(P(1)/(Kd+P(1))) - kdd*P(1);    P(3)*VM3*((1-P(2))/(K3+(1-P(2)))) - V4*(P(2)/(K4+P(2)));    VM1*(P(1)/(Kc+P(1)))*((1-P(3))/(K1+(1-P(3)))) - V2*(P(3)/(K2+P(3)))];
[T P] = ode45(PhsKnsKnt, [0  100],  [0.01;  0.01;  0.01]);
C = P(:,1);
X = P(:,2);
M = P(:,3);
figure(1)
plot(T, P)
legend('[C] (µ\itM\rm)',  '[X] (µ\itM\rm)',  '[M] (µ\itM\rm)', 'Location','NorthWest')
xlabel('Time (s)')
ylabel('Concentration')
grid
Just out of curiosity, what biochemical system is this?
Related Question