Need help converting berkeley madonna code to matlab code

Not sure where to start and how to convert the berkeley code to matlab

Here is the Berkeley Madonna code:

{Top model}

METHOD RK4

STARTTIME = 0

STOPTIME = 24

DT = 0.02

{Reservoirs Blood, Fat, NonFat and Liver}

d/dt (NF) = + Jnf

INIT NF = 0

d/dt (F) = + Jf

INIT F = 0

d/dt (B) = - Jnf - Jf - Jl + Jresp

INIT B = 0

d/dt (L) = + Jl - Jmetab

INIT L = 0

{Flows}

Jnf = Qnf*(Cb-Cnf)

Jf = Qf*(Cb-Cfv)

Jl = Ql*(Cb - Clv)

Jmetab = vmax*Cl/(Km + Cl)

{Replace squarepulse with " IF t>6h Jresp=0 "}

Jresp = Qp*(Ci - (Cb/Pb))*squarepulse(0,6)

{Functions}

Vnf = 1

Cnf = NF/Vnf

Qnf = 1

Vf = 1

Cf = F/Vf

Qf = 1

Cb = B/Vb

Vb = 1

Cfv = Cf/Pf

Pf = 20

Vl = 1

Cl = L/Vl

Ql = 1

vmax = 1

Km = 1

Pl = 2

Clv = Cl/Pl

Pb = 18

{Benzene conc in inhaled air}

Ci = 0.32

{alveolar ventilation rate}

Qp = 5.74

[Walter Roberson]

Perhaps you would want to code something like that in SimBiology ?

[Cris LaPierre]

Berkley Madonna is a differential equation solver. You can likely implement your code using MATLAB's main ode solver, ode45.

You might find these answers helpful.

[Walter Roberson]

SimBiology produces and solves differential equations -- it is one of the several options, one that might look closer to the Madonna code.

ode45() and related functions are other useful options for numeric solving. The ode functions whose name include 's' are useful for stiff systems.

It can also be useful to set up differential equations using the Symbolic Toolbox, and see if dsolve() can happen to solve them exactly, and if not then use odeFunction() to convert into functions for use with the numeric solvers.

Answer 1

[Cris LaPierre]

Here's equivalent MATLAB code.

STARTTIME = 0;
STOPTIME = 24;

% initial values
NF = 0;
F = 0;
B = 0;
L = 0;
y0 = [NF; F; B; L];

% solve ode
[t,y] = ode45(@odefun,[STARTTIME STOPTIME],y0);

% Plot results
yyaxis left
plot(t,y(:,[1,3]))
ylabel("NF,B")

yyaxis right
plot(t,y(:,[2,4]))
ylabel("F,L")

legend(["NF","B","F","L"])

function ddt = odefun(t,y)
NF = y(1);
F = y(2);
B = y(3);
L = y(4);

% Constants
Vnf = 1;
Cnf = NF/Vnf;
Qnf = 1;
Vf = 1;
Cf = F/Vf;
Qf = 1;
Vb = 1;
Cb = B/Vb;
Pf = 20;
Cfv = Cf/Pf;
Vl = 1;
Cl = L/Vl;
Ql = 1;
vmax = 1;
Km = 1;
Pl = 2;
Clv = Cl/Pl;
Pb = 18;
% Benzene conc in inhaled air
Ci = 0.32;
% alveolar ventilation rate
Qp = 5.74;

% Flows
Jnf = Qnf*(Cb-Cnf);
Jf = Qf*(Cb-Cfv);
Jl = Ql*(Cb - Clv);
Jmetab = vmax*Cl/(Km + Cl);

% Replace squarepulse with " IF t>6h Jresp=0 "
if t<=6
    Jresp = Qp*(Ci - (Cb/Pb));
else
    Jresp = 0;
end

% Reservoirs Blood, Fat, NonFat and Liver
ddt(1,1) = Jnf;
ddt(2,1) = Jf;
ddt(3,1) = -Jnf - Jf - Jl + Jresp;
ddt(4,1) =  Jl - Jmetab;

end

Compare those figures to the plots generated by BM

[Walter Roberson]

if t<=6
    Jresp = Qp*(Ci - (Cb/Pb));
else
    Jresp = 0;
end

That introduces a singularity. You need to stop the ode at t = 6 and restart it.

[Cris LaPierre]

True. There is no guarantee the solver will land anywhere near 6 otherwise.

[Arthur Goldsipe]

Oh, that might be one reason to check out SimBiology. It automatically restarts the solver after discontinuities, and it has a couple of different ways of implementing these discontinuities that might be more natural than these "if/else" blocks of code.

[Jeremy Huard]

SimBiology understand the concept of events to deal with those things without having to manually 'stop' the solver.

So, you could define your model in SimBiology so that:

• it contains an algebraic equation Jresp = Jrespon*Qp*(Ci - (Cb/Pb))
• Jrespon=1 at simulation start
• and an event with a trigger time>=6 and and event function Jrespon=0

Also, SimBiology could be useful if you want to do parameter estimation, use units, run sensitvitiy analyses, run parameter scans, etc.

Arthur pointed to a converter BM->SimBiology in his comment below if you would like to try SimBiology out.

And there is a series of video tutorials you might find useful:

SimBiology Tutorials for QSP, PBPK, and PK/PD Modeling and Analysis

[Jeremy Huard]

I'm attaching the SimBiology implementation of your model.

You can either open justing.sbproj in SimBiology and use the Model Analyzer App to run simulations (the file contains a pre-configured program to replicate Cris' plot), or you can use SimBiology commands to write scripts.

For instance:

% load model
sbioloadproject justin.sbproj
getequations(m1)

ans = 
    'ODEs:
     d(NF)/dt = Reaction_4
     d(F)/dt = Reaction_5
     d(B)/dt = Reaction_2 + Reaction_3 - Reaction_5 - Reaction_6
     d(L)/dt = -Reaction_1 + Reaction_6
     
     Fluxes:
     Reaction_1 = jmetab
     Reaction_2 = jresp
     Reaction_3 = -jnf
     Reaction_4 = jnf
     Reaction_5 = jf
     Reaction_6 = jl
     
     Repeated Assignments:
     cnf = NF/vnf
     cl = L/vl
     jmetab = vmax*cl/(km+cl)
     clv = cl/pl
     cf = F/vf
     cfv = cf/pf
     cb = B/vb
     jresp = qp*(ci-(cb/pb))*jresp_on
     jnf = qnf*(cb-cnf)
     jl = ql*(cb-clv)
     jf = qf*(cb-cfv)
     
     Parameter Values:
     ci = 0.32
     km = 1
     pb = 18
     pf = 20
     pl = 2
     qf = 1
     ql = 1
     qnf = 1
     qp = 5.74
     vb = 1
     vf = 1
     vl = 1
     vmax = 1
     vnf = 1
     System = 1
     cb = 0
     cf = 0
     cfv = 0
     cl = 0
     clv = 0
     cnf = 0
     jf = 0
     jl = 0
     jmetab = 0
     jnf = 0
     jresp = 1.8368
     
     Initial Conditions:
     NF = 0
     F = 0
     B = 0
     L = 0
     jresp_on = 1
     
     Events:
     Name                          	Trigger                       	Function                      
                                   	time >= 6                     	jresp_on = 0                  
     '

% run simulation
sd = sbiosimulate(m1);
[t,y] = selectbyname(sd, ["NF","B","F","L"]);

% Plot results
figure;
ax = axes(XLimitMethod='padded', YLimitMethod='padded');
yyaxis(ax, 'left');
plot(t,y(:,[1,2]))
ylabel("NF,B")
yyaxis(ax, 'right');
plot(t,y(:,[3,4]))
ylabel("F,L")
grid(ax,'on');
legend(["NF","B","F","L"],Location='eastoutside',Box='off')

I hope this helps.

[Walter Roberson]

@Jeremy Huard

Those sound like useful features of SimBiology.

In Answers, we get a fair number of questions about ode45() and ode23s() in which people use if/else in discontinuous ways, and we have to keep talking to them about how that is (usually) mathematically incompatible with the ode*() functions. Is SimBiology an alternative we should be talking up for general ODE systems that have events and impulses?

... and would there just happen to be a conversion function that could take symbolic ODE with dirac or heaviside or piecewise and build and run an appropriate SimBiology system?

Answer 2

[Arthur Goldsipe]

If you'd like to convert a Berkeley Madonna model to the equivalent SimBiology model, you can try using this converter from the File Exchange.