Need help converting berkeley madonna code to matlab code

Question

Justin Lee 2021-5-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/833543-need-help-converting-berkeley-madonna-code-to-matlab-code

评论： Walter Roberson 2021-5-20

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Cris LaPierre 2021-5-18

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 2021-5-18

编辑：Walter Roberson 2021-5-18

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 2021-5-18

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/833543-need-help-converting-berkeley-madonna-code-to-matlab-code#answer_702823

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Jeremy Huard 2021-5-19

justin.sbproj

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 2021-5-20

@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?

请先登录，再进行评论。