Simbiology with compartments of different volume

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)

Answer 1

[Arthur Goldsipe]

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.

[felix MP]

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]

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.