Trying to write a 3 compartment model
7 views (last 30 days)
Show older comments
% set parameter values
Vd = 210; % volume of distribution litres = vol compartment 2
V1 = 10; % volume of gut compartment
V3 = 30;
Dose = 25.*1000; % milligrams* 1000 gives micrograms
C0 = Dose/V1;
k1 = 0.24; % per hour
k2 = 0.28; % per hour
vmax = .47;
km = .38;
a = .42;
tmax = 12; % time range over which to solve model (hours)
Cinit = [C0,0,0]; % initial value for compartment 1 and compartment 2
% write down the expressions for the fluxes
J1i = @(t,C) 0; % no continuous influx to compartment 1
J1e = @(t,C) k1.*V1.*C(1); % efflux from compartment 1
J2i = J1e;
J2e = @(t,C) k2.*V2.*C(2); % influx to compartment 2 = efflux from compartment 1
J3i = J2e;
J3e = @(t,C) k2.*V3.*C(3); % efflux from compartment 3
% solve the differential equation
sol = ode45(@ThreeCompartmentModel,[0 tmax],Cinit,[],J1i,J1e,J2i,J2e,J3i,J3e,V1,Vd,V3);% plot results
% plot results
tstep = 0.01;
tspan = 0:tstep:tmax;
y = deval(sol,tspan);
plot(tspan,y)
legend('Compartment 1','Compartment 2','Concentration 3')
ylabel('drug concentration')
xlabel('time (hours)')
I wanted to do this in SimBiology or something. And be able to do multiple doses. But when I look at the SimBiology, it is modeling data I think. Can I simply represent a system like this? I wanted to model this. https://php.radford.edu/~ejmt/deliveryBoy.php?paper=eJMT_v1n3p3
2 Comments
Praveen Kumar M
on 19 Nov 2018
Hi Ryan,
I am clinical pharamcologist. Have you found the answer? If yes, please share the solution with us. And if you have not found the answer, then we can combinely work towards the solution. You can contact me either through MathWorks or by gmail (my mail: praveenkumarpgiindia@gmail.com)
Praveen Kumar M
on 19 Nov 2018
In the above code, which you have shared, I have a query. How have you defined the ThreeCompartmentModel function or variable?
Answers (2)
Achal Mahajan
on 1 Jun 2022
Hi Ryan,
Here is the implementation of the problem in SimBiology. The following implementation is motivated from the paper that you have shared. To understand it better, as Arthur have mentioned above, please use the examples/videos/webinars.
clear all;clc;close all;
model = sbiomodel('Q4_MATLAB');
comp1Obj = addcompartment(model,'Stomach');
comp1Obj.CapacityUnits = 'liter';
comp1Obj.Value = 1;
comp2Obj = addcompartment(model,'Small_Int');
comp2Obj.CapacityUnits = 'liter';
comp2Obj.Value = 1;
comp3Obj = addcompartment(model,'Central');
comp3Obj.CapacityUnits = 'liter';
comp3Obj.Value = 1;
% Let's define all the reactions in their respective compartments
r1 = addreaction(model,'Stomach.Alcohol -> Small_Int.Alcohol');
r1.ReactionRate = 'k1*Stomach.Alcohol/(1+a*Stomach.Alcohol^2)';
r2 = addreaction(model,'Small_Int.Alcohol -> Central.Alcohol');
k = r2.addkineticlaw('MassAction');
k.ParameterVariableNames = 'k2';
r3 = addreaction(model,'Central.Alcohol -> null');
k = r3.addkineticlaw('Henri-Michaelis-Menten');
k.ParameterVariableNames = {'Vmax','km'};
k.SpeciesVariableNames = 'Central.Alcohol';
% Initial concentration and specify units
for i = 1:length(model.Species)
model.Species(i).Units = 'gram/liter';
end
model.Species(1).InitialAmount = 0.455;% This is what is mentioned in the paper, different from what you have implemented
model.Species(2).InitialAmount = 0;
model.Species(3).InitialAmount = 0;
% Now we can add the parameters of the model like rate constants
addparameter(model,'Vmax','Value',0.470,'ValueUnits','gram/liter/hour');
addparameter(model,'km','Value',0.380,'ValueUnits','gram/liter');
addparameter(model,'k1','Value',5.55,'ValueUnits','1/hour');
addparameter(model,'k2','Value',7.05,'ValueUnits','1/hour');
addparameter(model,'a','Value',0.42,'ValueUnits','liter^2/gram^2');
% Let's run the simulation now
cs = getconfigset(model,'active');
cs.StopTime = 4;
cs.TimeUnits = 'hour';
% Run the simulation and plot the data
cs.CompileOptions.UnitConversion = true;
[t, simdata, names] = sbiosimulate(model);
%[t, simdata, names] = sbiosimulate(model,d1); %Uncomment this if you want to add a dose
%[t, simdata, names] = sbiosimulate(model,v2); % or v3 Uncomment this if you want to add a variant
plot(t,simdata)
legend(names,'Location','NorthEastOutside');
-Achal
0 Comments
See Also
Categories
Find more on Scan Parameter Ranges in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!