Simbiology with compartments of different volume
7 次查看(过去 30 天)
显示 更早的评论
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)
0 个评论
回答(1 个)
Arthur Goldsipe
2023-2-21
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 个评论
Arthur Goldsipe
2023-2-22
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.
社区
更多回答在 SimBiology Community
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Scan Parameter Ranges 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!