Implementing a code from Berkley Madonna into Matlab

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.

[Star Strider]

First: No. It’s not even close. If you post the equations you want to integrate, and the code you wrote to integrate them, we can help you if you have problems.

Second, if you simply wrote your own functions with the same names as the MATLAB functions, all may not be lost. (This is known as ‘shadowing’ MATLAB functions.) See the documentation for the which function for details on how to determine this.

[Ingrid Tigges]

Are you aware of this Import Berkley Madonna models into SimBiology entry from the file exchange? It gives you an easy way to implement the models into SimBiology.

[Star Strider]

That’s good to know. Thank you.

KayLynn doesn’t have SimBiology, though. (That was in a subsequent post.) Neither do I.

Answer 1

[Star Strider]

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?

[KayLynn]

Thanks for the start it makes some sense as to how you laid it out. The only information I was given was that it is a Mitotic Oscillator.

[Star Strider]

My pleasure!

It definitely displays periodic behaviour. I experimented with both ode45 and ode15s. The ode15s solver is probably more appropriate for it.

I’ll look up ‘Mitotic Oscillator’. That’s new. I’m more into intermediary metabolism, so I didn’t recognise it.

[Star Strider]

I searched on ‘Mitotic Oscillator’, and came across a page documenting the MATLAB SimBiology implementation of it. (I don’t have SimBiology.) If you have access to SimBiology, you might want to give this page a look: Minimal Cascade Model for a Mitotic Oscillator. It also has a thorough discussion of the model.

Thanks also for the follow-up. It turns out I actually read that paper a few years after it was published, but I didn’t recognise the model from the equations.

[KayLynn]

are you familiar with perturbations and dtermining SFO rankings??

[KayLynn]

Hi star strider....since you have been so helpful to me on my question yesterday...I have one more if you wouldnt mind providing some insight...it deals with the paper I mentioned yesterday for the Mitotic Oscillator..... I now want to "break the system" to disrupt the oscillations (mainly looking at the discussion section of the paper)....I may be missing the concept completely but Im not sure how to do this in matlab code without completely screwing the entire code...any help you could provide would be appreciated. THanks

[Star Strider]

I missed your earlier post about perturbations and SFO rankings. I’m not familiar with them, but they look like numerical derivatives. I’ll do a search later to see what I can find. (I found it interesting that when I plotted the ratio of the two for Kd, the effect was most pronounced at lower concentrations. If Kd is the dissociation constant, that entirely makes sense.)

I looked at both threads and didn’t find an uploaded paper or a mention of it. (The only reference I saw was my link to the SimBiology documentation on the mitotic oscillator. It linked to two PNAS papers by Goldbeter, and I believe are freely available online. Did you mean one of them?) If you upload the PDF of the paper you referred to as an attachment to your initial post here, I’ll be glad to take a look at it.

You probably want to ‘poison’ the system somehow, but I’m not having read the paper yet, I’m not certain how you want to do that.

I’ll look for the uploaded paper (or a link to the free paper PDF) and reply when I’ve read it.

[KayLynn]

thank you...i think it deals with maybe the time points but I am not com;etely certain with that

[Star Strider]

I need to be offline for a bit, but I’ll download it as soon as I see you’ve uploaded it. Not having read it, I can’t comment further just now.

[Star Strider]

I couldn’t find anything on ‘SFO ranking’ in an Internet search (DuckDuckGo and PubMed). The problem with ‘perturbation’ is that it means different things in different contexts. I couldn’t find anything dealing specifically with enzyme kinetics.

As far as ‘poisoning’ the system, that probably only requires changing one or more binding constants. Shouldn’t screw up the code for the DEs.

This paper ‘Arresting the mitotic oscillator and the control of cell proliferation: insights from a cascade model for cdc2 kinase activation’ may be what you want to do, but it’s not available online.