Simbiology: Generate reactions with a loop

I'm creating a pk model with the following reaction rate: 'Dose*ka*e^(-ka*(time-TimeOfDose))'. I haven't been able to figure out how to do this in Simbiology, other that to just create an absorbance reaction for each dose. However, I need a lot of doses, and to be able to easily change the number of doses. I figured that I would use a loop to generate the correct number of reactions, however I'm getting some errors when trying this.

Here is the code:

for i = 0:(nDose - 1)

   r(i) = addreaction(m1, 'null -> Central.Drug');
   set(r(i), 'ReactionRate', Dose*ka*e^(-ka*(time-(i*7)));
end

With 7 being the dosing interval. This gives the error:

 Subscript indices must either be real positive integers or logicals.

Answer 1

[Arthur Goldsipe]

The immediate problem is that MATLAB uses 1-based indexes when indexing into vectors, so you need to make sure that you start at r(1) instead of r(0).

But it sounds like you're trying to model first-order absorption of a drug. There is a simpler way to do that in SimBiology. Basically, add an intermediary species and reaction to your model to account for the first order absorption. For example, if you call the species Dose and put the species in compartment Central, then add a reaction "Central.Dose -> Central.Drug" with reaction rate "ka*Central". Finally, implement your dosing using a dose object with "Dose" as the target.

[Colin Cess]

Sorry, I've been pretty busy and just had time to look at this.

I don't think I was too clear in my question, but the equation that I have written is only part of the entire absorption equation that I'm using (I didn't write out the entire thing since the equation doesn't matter too much, just the concept). I'm modeling a subcutaneous dose, and what I wrote is the first absorption stage. There are two more stages, using transit compartments, and the diffeq for the central compartment is 'explicit solution for absorption - elimination'. I've been able to simulate a single dose just fine, but multi-dose is giving me some trouble. For each dose, I need to use the absorption equation where time t = time-TimeOfDose since I'm not adding new doses into the same compartment, but I haven't been able to figure out how to get this to work easily. I know I could just make a reaction for each dose, but that would be between 15 and 30 doses, and it wouldn't be easy to change the number of doses.

Hopefully this is more clear.

[Arthur Goldsipe]

I guess I still don't understand why you want to use the analytical/explicit solution to differential equations with SimBiology. I'll show you one way that uses events below. But I think it's going to be a lot more complicated than just letting SimBiology solve the equations for you. If you have additional follow-up questions, I suggest contacting me through my profile page. I don't get notified when you post a comment here.

m1 = sbiomodel('m1');
m1.addparameter('Dose', 2.0);
m1.addparameter('ka', 0.1);
m1.addparameter('doseTerm', 0.0, 'ConstantValue', false);
m1.addrule('doseTerm = Dose', 'initialAssignment'); % t0 dose
r1 = m1.addreaction('null -> Central.Drug'); 
r1.ReactionRate = 'doseTerm*ka*exp(-ka*time)';
m1.addparameter('doseInterval', 7);
m1.addparameter('doseCounter', 1, 'ConstantValue', false);
m1.addparameter('nDose', 3);
m1.addevent('time > doseCounter*doseInterval && doseCounter < nDose', ...
    {'doseTerm = doseTerm + Dose*ka*exp(ka*doseCounter*time)', ...
    'doseCounter = doseCounter + 1'});