# Need help converting berkeley madonna code to matlab code

16 views (last 30 days)
Justin Lee on 18 May 2021
Commented: Walter Roberson on 20 May 2021
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 on 18 May 2021
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.

Cris LaPierre on 18 May 2021
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 on 20 May 2021
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?

Arthur Goldsipe on 18 May 2021
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.