Simbiology with compartments of different volume

7 views (last 30 days)
Hi there,
I stumbled upon a (for me ) non-intuitive simulation behaviuor. Consider two compartments of differnet volume, intra and extra, and one species S that can freely diffuse across the compartments. I set it up programmatically and tried to correct for the volume difference with a clearance parameter as described here https://nl.mathworks.com/help/simbio/ug/creating-pharmacokinetic-models.html
Initial conc are S = 0 in the intra comaprtment and S=1 in the extra compartment.
Results vary with the intra compartment volume, and seem to make sense unless the intra compartment gets very small (1E-10): What I found very weird is that the intra.S cocnentration then increased linearly (!) from 0 to 1.2 (!) with a slope that is directly proportional to the simulation end time (!).
2 Questions:
1) Is my clearance implementation correct at all or how do i do this better?
2) Why do i get these non intuitive results if one compartment becomes small?
Below is the script. Play with the simulation end time (simuTime) and intracellComp.Value
Thanks in advance! Felix
clc
clear
plotspeciesnames = {'extra.S', 'intra.S'};
InitialAmount.intra.S = 0; % units: M
InitialAmount.extra.S = 1; % units: M
clearance = 0.1; % units: volume / time
simuTime = 3000;
%% set up simbio model with 2 compartments, extra and intra
model = sbiomodel('txn');
extracellComp= addcompartment(model,"extra");
extracellComp.Value = 1;
extracellComp.Units= 'liter';
intracellComp= addcompartment(model,"intra");
intracellComp.Value = 1e-10;
intracellComp.Units= 'liter';
%% Influx of nucleotides
kin = clearance / extracellComp.Value;
kout = clearance / intracellComp.Value;
reactionText = 'extra.S <-> intra.S';
influx = addreaction(model, reactionText);
kineticLaw = addkineticlaw(influx, 'MassAction');
addparameter(model, 'kin', kin, 'ValueUnits','1/minute');
addparameter(model, 'kout', kout, 'ValueUnits','1/minute');
kineticLaw.ParameterVariableNames = {'kin', 'kout'};
%% Set units of all species to molarity
species=get(model, 'Species');
set(species, 'Units', 'molarity')
%% Set up initial conditions
species = sbioselect(intracellComp, 'Name', 'S');
species.InitialAmount = InitialAmount.intra.S;
species = sbioselect(extracellComp, 'Name', 'S');
species.InitialAmount = InitialAmount.extra.S;
%% display all reactions and reaction rates
get(model.Reactions, {'Reaction', 'ReactionRate'})
get(model, 'Parameters')
get(model, 'Species')
%% Set up solver
configset = model.getconfigset;
configset.CompileOptions.UnitConversion = true;
configset.StopTime = simuTime;
configset.CompileOptions.DimensionalAnalysis = true;
%% Simulate the model
sim = sbiosimulate(model);
%% plot results
for ii=length(plotspeciesnames):-1:1
simdataSingleSpecies = selectbyname(sim, plotspeciesnames{ii} );
simdata(:,ii) = simdataSingleSpecies.Data;
end
time = simdataSingleSpecies.Time;
plot(time, simdata)
legend(plotspeciesnames)

Answers (1)

Arthur Goldsipe
Arthur Goldsipe on 21 Feb 2023
I think this is because of the volume scaling of the reaction rate. I suggest calling getequations on your model, and seeing if the resulting reaction rate expression is what you intended. You should see the following in the output:
Reaction_1 = (kin*extra.S)*extra-(kout*intra.S)*intra
Note that the forward reaction is scaled by the volume extra, while the reverse reaction is scalled by the volume intra. This sort of scaling only happens in SimBiology if you select reversible mass action kinetics. If this is not what you want, then one option would be to set the kinetic law to "Unknown" and set the reaction rate to the exact expression you want.
SimBiology needs to convert concentrations to amounts before solving the equations. When your reaction rate has dimensions of concentration per time, then SimBiology will multiple the reaction rate by an appropriate volume. Under most circumstances, only one volume gets used for this scaling. But in the special case of reversible mass action kinetics, SimBiology scaled the forward reaction rate to be scaled by the reactant's compartment, and the reverse reaction rate to be scaled by the product's compartment. Such scaling is consistent with the underlying assumptions of mass action kinetics.
  2 Comments
felix MP
felix MP on 21 Feb 2023
Thanks, Arthur. Not sure if that's the problem though (as i tried to correct for it already using rate conmstants that scale inversely with the compartment volumes).
Rather, in my own testing the massiv (!) deviations from expected results appear to be a sensitivity issue of the solver. Try running the code below with or without the line
configset.SolverOptions.MaxStep = 0.1;
Bewildering. I guess that begs the question how i can ever be sure that simulated results are reliable. Any best preactices?
Felix
With the MaxStep line:
clc
clear
plotspeciesnames = {'extra.S', 'intra.S'};
InitialAmount.intra.S = 0; % units: M
InitialAmount.extra.S = 1; % units: M
clearance = 0.1; % units: volume / time
simuTime = 3000;
%% set up simbio model with 2 compartments, extra and intra
model = sbiomodel('txn');
extracellComp= addcompartment(model,"extra");
extracellComp.Value = 1;
extracellComp.Units= 'liter';
intracellComp= addcompartment(model,"intra");
intracellComp.Value = 1e-10;
intracellComp.Units= 'liter';
%% Influx of nucleotides
kin = clearance / extracellComp.Value;
kout = clearance / intracellComp.Value;
reactionText = 'extra.S <-> intra.S';
influx = addreaction(model, reactionText);
kineticLaw = addkineticlaw(influx, 'MassAction');
addparameter(model, 'kin', kin, 'ValueUnits','1/minute');
addparameter(model, 'kout', kout, 'ValueUnits','1/minute');
kineticLaw.ParameterVariableNames = {'kin', 'kout'};
%% Set units of all species to molarity
species=get(model, 'Species');
set(species, 'Units', 'molarity')
%% Set up initial conditions
species = sbioselect(intracellComp, 'Name', 'S');
species.InitialAmount = InitialAmount.intra.S;
species = sbioselect(extracellComp, 'Name', 'S');
species.InitialAmount = InitialAmount.extra.S;
%% display all reactions and reaction rates
get(model.Reactions, {'Reaction', 'ReactionRate'})
ans = 1×2 cell array
{'extra.S <-> intra.S'} {'kin*extra.S - kout*intra.S'}
get(model, 'Parameters')
ans =
SimBiology Parameter Array Index: Name: Value: Units: 1 kin 0.1 1/minute 2 kout 1e+09 1/minute
get(model, 'Species')
ans =
SimBiology Species Array Index: Compartment: Name: Value: Units: 1 extra S 1 molarity 2 intra S 0 molarity
getequations(model)
ans =
'ODEs: d(extra.S)/dt = 1/extra*(-Reaction_1) d(intra.S)/dt = 1/intra*(Reaction_1) Fluxes: Reaction_1 = (kin*extra.S)*extra-(kout*intra.S)*intra Parameter Values: kin = 0.1 kout = 1e+09 extra = 1 intra = 1e-10 Initial Conditions: extra.S = 1 intra.S = 0 '
%% Set up solver
configset = model.getconfigset;
configset.SolverOptions.MaxStep = 0.1;
configset.CompileOptions.UnitConversion = true;
configset.StopTime = simuTime;
configset.CompileOptions.DimensionalAnalysis = true;
%% Simulate the model
sim = sbiosimulate(model);
%% plot results
for ii=length(plotspeciesnames):-1:1
simdataSingleSpecies = selectbyname(sim, plotspeciesnames{ii} );
simdata(:,ii) = simdataSingleSpecies.Data;
end
time = simdataSingleSpecies.Time;
plot(time, simdata)
legend(plotspeciesnames)
Without the MaxStep line:
clc
clear
plotspeciesnames = {'extra.S', 'intra.S'};
InitialAmount.intra.S = 0; % units: M
InitialAmount.extra.S = 1; % units: M
clearance = 0.1; % units: volume / time
simuTime = 3000;
%% set up simbio model with 2 compartments, extra and intra
model = sbiomodel('txn');
extracellComp= addcompartment(model,"extra");
extracellComp.Value = 1;
extracellComp.Units= 'liter';
intracellComp= addcompartment(model,"intra");
intracellComp.Value = 1e-10;
intracellComp.Units= 'liter';
%% Influx of nucleotides
kin = clearance / extracellComp.Value;
kout = clearance / intracellComp.Value;
reactionText = 'extra.S <-> intra.S';
influx = addreaction(model, reactionText);
kineticLaw = addkineticlaw(influx, 'MassAction');
addparameter(model, 'kin', kin, 'ValueUnits','1/minute');
addparameter(model, 'kout', kout, 'ValueUnits','1/minute');
kineticLaw.ParameterVariableNames = {'kin', 'kout'};
%% Set units of all species to molarity
species=get(model, 'Species');
set(species, 'Units', 'molarity')
%% Set up initial conditions
species = sbioselect(intracellComp, 'Name', 'S');
species.InitialAmount = InitialAmount.intra.S;
species = sbioselect(extracellComp, 'Name', 'S');
species.InitialAmount = InitialAmount.extra.S;
%% display all reactions and reaction rates
get(model.Reactions, {'Reaction', 'ReactionRate'})
ans = 1×2 cell array
{'extra.S <-> intra.S'} {'kin*extra.S - kout*intra.S'}
get(model, 'Parameters')
ans =
SimBiology Parameter Array Index: Name: Value: Units: 1 kin 0.1 1/minute 2 kout 1e+09 1/minute
get(model, 'Species')
ans =
SimBiology Species Array Index: Compartment: Name: Value: Units: 1 extra S 1 molarity 2 intra S 0 molarity
getequations(model)
ans =
'ODEs: d(extra.S)/dt = 1/extra*(-Reaction_1) d(intra.S)/dt = 1/intra*(Reaction_1) Fluxes: Reaction_1 = (kin*extra.S)*extra-(kout*intra.S)*intra Parameter Values: kin = 0.1 kout = 1e+09 extra = 1 intra = 1e-10 Initial Conditions: extra.S = 1 intra.S = 0 '
%% Set up solver
configset = model.getconfigset;
%configset.SolverOptions.MaxStep = 0.1;
configset.CompileOptions.UnitConversion = true;
configset.StopTime = simuTime;
configset.CompileOptions.DimensionalAnalysis = true;
%% Simulate the model
sim = sbiosimulate(model);
%% plot results
for ii=length(plotspeciesnames):-1:1
simdataSingleSpecies = selectbyname(sim, plotspeciesnames{ii} );
simdata(:,ii) = simdataSingleSpecies.Data;
end
time = simdataSingleSpecies.Time;
plot(time, simdata)
legend(plotspeciesnames)
Arthur Goldsipe
Arthur Goldsipe on 22 Feb 2023
Ah, it looks like I jumped to conclusions too quickly. Also, thanks for confirming the equations are what you expect.
When you have such extreme differences in compartment volumes, the differential equations become very stiff (making them hard to solve). You need to be extra careful in such situations and inspect your results (just like you're doing). In such situations, the absolute tolerance scaling algorithm can break down. You might want to disable that option and see if your results are more reasaonable.
You can read more about tolerances here, and you can find some simulation troubleshooting tips here.
As for best practices, the first thing you should always do is see if your model inputs make sense. Is it reasonable to expect the compartments in your model to differ by a factor of 10^10? Next, make sure your model outputs make sense. In cases like this where they don't, then you need to look deeper. If you don't know whether the outputs make sense and you want to make sure there's not a problem, you can do further analysis. You might want to see whether the results are extremely sensitive to your initial conditions or your current solver conditions. If they are extremely sensitive, then it may indicate that your model is extremely difficult to solve accurately. If that's the case, then it might make sense to see whether you should change your model. In the example you provide where two compartments have very different volumes, perhaps it makes sense to treat the concentration in the larger compartment as effectively constant, and only model the concentration change in the smaller compartment.

Sign in to comment.

Communities

More Answers in the  SimBiology Community

Categories

Find more on Extend Modeling Environment in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!