Simbiology dose an ode

Hello,

I have a question concerning of the use of simbiology adddose (or sbiodose) commands. I have defined a compartment, and four species (that will be my state variables in the model). I have also defined a bunch of parameters and added to the model. After this step, I have explicitly defined four differential equations that describe my model. (I did it in this way since I have no proper background in reactons and stuff). My goal was to simulate the system with a prespecified dosing scheulde that came from an experiment. When i create my dose and want to specify the target name as 'x3', matlab throws an error: 'Invalid dose target 'x3' in dose 'd'. This object cannot be the Rule variable of any rule'. I checked the documentation and forums but I haven't found anything related to my issue. Any help much appreciated, code supplemented below.

[a1,b1,c1,ED501,k11,k21,n1,w1,x10] = getParams(1); %This get some exact values

%This part of the code gets the injection times and amounts from the experiment
load('measurements.mat');
measurement = measurement{1};
mouseID = measurement.vMouseID{1};
tInjections = measurement.tMeasurement{1};
vMeasuredVolumes = measurement.vVolumePhysCon{1};
vInjectionDose = measurement.vDose{1};

doseIdx = find(vInjectionDose ~= 0);

m = sbiomodel('m');
comp = addcompartment(m,'comp');

%State variables, input, output
x1 = addspecies(m,'x1','InitialAmount',x10); %Living
x2 = addspecies(m,'x2','InitialAmount',0); %Dead
x3 = addspecies(m,'x3','InitialAmount',0); %Conc
x4 = addspecies(m,'x4','InitialAmount',0); %Periph Conc
u = addspecies(m,'u','InitialAmount',0);
y = addspecies(m,'y','InitialAmount',x10);

%Model parameters
a = addparameter(m,'a','Value',a1);
b = addparameter(m,'b','Value',b1);
n = addparameter(m,'n','Value',n1);
w = addparameter(m,'w','Value',w1);
ED50 = addparameter(m,'ED50','Value',ED501);
c = addparameter(m,'c','Value',c1);
k1 = addparameter(m,'k1','Value',k11);
k2 = addparameter(m,'k2','Value',k21);

%Differential equations
dxdt1 = addrule(m,'x1 = (a-n) * x1 - b * ((x1*x3)/(ED50 + x3))','RuleType','rate');
dxdt2 = addrule(m,'x2 = n * x1 + b*((x1*x3)/(ED50 + x3)) - w*x2','RuleType','rate');
dxdt3 = addrule(m,'x3 = -(c + k1)*x3 + k2*x4','RuleType','rate');
dxdt4 = addrule(m,'x4 = k1 * x3 - k2*x4','RuleType','rate');
y = addrule(m,'y = x1 + x2','RuleType','repeatedAssignment');

%Add dose - this is where things go south for me
d = sbiodose('d','schedule');
d.Amount = vInjectionDose(doseIdx);
d.Time = tInjections(doseIdx);
d.TargetName = 'comp.x3';
d.Active = 1;
[t,x,names] = sbiosimulate(m,d);

plot(t,x(:,1:4))

Best regards,

Bence

EDIT: I forgot to mention, that my goal was to direcly inject the drug into species 'x3' with an impulsive action (bolus administration), but the ode solver does not permit such, since 'x3' is a rate variable. My guess is that this can be possible somehow by defining it as a reaction, but dxdt3 and dxdt4 is the clearence equations that I want x3 to satisfy.

Answer 1

[Jeremy Huard]

Hi Bence,

you're right, doses can't be applied to species which are defined by rate rules. Their dynamics must be defined by reactions instead.

In your specific case, converting the rate rules into reactions is pretty quick. You can replace all dxdt definitions by the following:

addreaction(m, 'null -> x1', 'ReactionRate','a * x1');
addreaction(m, 'x1 -> x2', 'ReactionRate','n * x1');
addreaction(m, 'x1 + x3 -> x2 + x3', 'ReactionRate','b * ((x1*x3)/(ED50 + x3))');
addreaction(m, 'x2 -> null', 'ReactionRate','w*x2');
addreaction(m, 'x3 <-> x4', 'ReactionRate','k1*x3 - k2*x4');
addreaction(m, 'x3 -> null', 'ReactionRate','c*x3');

If you don't use units, SimBiology will assume ba default that all species have dimensions of concentrations. That means that the ODEs derived from the reactions above will have a term (1/comp) in the right-hand side of the equations. To avoid this you can tell SimBiology to consider the species as amounts instead of concentrations.

You can insert the following lines before your call to sbiosimulate. Nothing else needs to be changed.

cs = getconfigset(m);
cs.CompileOptions.DefaultSpeciesDimension = 'substance';

Best regards,

Jérémy

[Ben]

Hi,

Thank you Jérémy, it works like a charm!

Best,

Bence

[Jeremy Huard]

One quick comment:

You can visualize the ODEs generated from the reactions with:

getequations(m)

You can also import your model into the SimBiology apps to visualize / edit/ analyze your model:

simbiologyModelBuilder(m)

[Jeremy Huard]

Oh, last thing:

I just realized that your defined your dose as Active.

This means that the dose will always be applied even if you don't pass it to sbiosimulate.

I suggest to leave the Active property to false and only apply the dose consciously when needed.

[Ben]

Thank you!