Time-dependent Growth Rate of Species in Model

Hello,

I'm trying to implement in Simbiology the model of a system which describes the kinetics of many species. The kinetics of one these species is non-standard in that it's change rate is time dependent with different kinetics in three different periods. In the first period, which occurs in the first five days, the species is able to distribute from its initial compartment to another at a fixed rate. Following the distribution of the species into the second compartment, in the second period which occurs for the next seven days, the species grows at a fixed constant rate. Finally, for the next seven days during the third period, the species decays at a fixed constant rate. Thus, if I'm not mistaken, the species has a time-dependent growth rate that can be described as a step function in time with changes in time = 5, 12, and 19 days. Is there a way implement these kinetics using Simbiology?

Many thanks in advance for any suggestions or comments.

Answer 1

[Florian Augustin]

Hi Marco,

you may be able to use SimBiology Events to switch between reaction rates. The general idea would be as follows:

• Define a reation rate rate1, rate2, rate3, ... for each of the phases you describe.
• Add dummy (phase-indicator) parameters phase1, phase2, phase3, ... to your model. The initial value of phase1 should be 1, all other parameters should have an initial value of 0.
• Add SimBiology Events to switch between different phases.

Below is an example using SimBiology on the MATLAB Command Window. Dependent on your actual rates you may be able to condense/rewrite the events; you may not even need the phase-indicator parameters. But I hope this demonstrates how you could achieve different reation rates in different time intervals.

% Define model

model = sbiomodel("time-dependent reaction rates");

addspecies(model, "A", 10);

addparameter(model, "k", 1, "Constant", false);

% Add phase-indicator (dummy) parameters:

addparameter(model, "phase1", 1, "Constant", false);

addparameter(model, "phase2", 0, "Constant", false);

addparameter(model, "phase3", 0, "Constant", false);

% Define reaction rate: here I am using rate1 = k*A, rate2 = -k, rate3 = -k*A as an example.

reaction = addreaction(model, "A -> null");

reaction.ReactionRate = "phase1*k*A - phase2*k - phase3*k*A";

% Add events to switch between phases:

addevent(model, "time >= 5" , ["phase1 = 0", "phase2 = 1", "k = 0.1"]);

addevent(model, "time >= 12", ["phase2 = 0", "phase3 = 1", "k = 0.4"]);

% Simulate model

configset = getconfigset(model);

configset.StopTime = 19;

configset.RuntimeOptions.StatesToLog = "A";

simData = sbiosimulate(model);

% Plot results

sbioplot(simData)

Best,

Florian

[Florian Augustin]

Instead of using a reaction, you could also use a SimBiology Rate Rule. This may make changing rates more "readable" when looking at the SimBiology model (in case the kinetic is not as simple as it is in my example).

[mitpi_182]

Hello Florian,

Many thanks for your response. It's very clearly written and answers my question effectively. One quick follow-up question: when one defines the event functions, such as '"time > 5", how does one define the units of time? Specifcally, I want to make sure that "time > 5" means time greater than five days. Also, I'm using the Simbiology GUI, is there anything that needs to done differently in there? Thanks again for your time!

[Florian Augustin]

Hi Marco,

The time units are specified in the Simulation Settings. You can access those from the toolstrips of the SimBiology Model Builder and SimBiology Model Analyzer apps. For example:

In the SimBiology Analyzer app you can also override those time units in individual programs. For example:

You can do unit conversion either manually, or let Simbiology do it for you. Assuming TimeUnits is specified as hour and you want the event to trigger after 5 days, then you would have to specify the event trigger as time >= 120. To let SimBiology do the unit conversion for you, you have to turn on UnitConversion (see Simulations Settings). You can then use parameters with specified units in your event triggers:

1. Add a parameter START_PERIOD with value 5 and unit day to your model. then you can write your event function as follows: time >= START_PERIOD.
2. Add a "unit" parameter [day] with value 1 and unit day to your model, then you can write your event function as follows: time >= 5*[day]. The square brackets have no meaning here. I just like this notation because it identifies this parameter as a "unit" converter, rater than a model parameter.

I hope this helps.

Best,

Florian

[mitpi_182]

Hello,

Hopefully its okay that I revive this thread. I recently tried to implement the kinetics I described. I followed the instructions discussed here to the dot. However, for some reason I can't get the simulation to run. I keep getting the following error: "Invalid input argument of type double."Input must be a structure or a Java or COM object." I was wondering if there are any common mistakes that trigger this error. I've looked everywhere for the bug and have not found anything. Any comments or suggestions are greatly welcome.

[Florian Augustin]

Hi Marco,

One problem could be that I said to add a parameter [day]. Can you name the parameter to day (without the square brackets) and see if the error persists? You can still refer to the parameter in events/expressions as [day], no changes needed there.

If you continue to see this error, would you be able to share the sbproj / script file you are working with? If you cannot share your model publicly on MATLAB Answers, please contact MathWorks Technical Support (please reference this post here in the bug report).

Thank you,

Florian

[mitpi_182]

Hello Florian,

Many thanks for your response. I didn't include the parameter expression you mentioned but still got the error. I decided to start re-implement the kinetics from scratch and was able to simulate the kinetics I needed. I'm not sure why I did not get the error this time. It'd be good to know for future reference so I'll send the file to technical support for analysis.

I have one final question to ask regarding this model. Hopefully it is okay that I ask this here and not at another thread. As described earlier, the kinetics that I'm seeking display two sharp transitions: one where proliferation occurs starting from day 5 and then one where degradation occurs seven days after that. I've run some simulations and believe I'm getting the right result. If plotted on a semi-log scale, the T cell kinetics look correct: at the onset of proliferation we get a straight line with positive slope followed by a straight line with negative slope at the onset of degradation. However, the transitions are quite sharp rendering the profiles to be non-differentiable at the two transition points. This is not how the profile looks in the figure of the example you provided, despite having similar sharp transitions in the kinetics. Is was wondering if there are any techniques we could use to smooth things out at the transition points using Simbiology?

Many thanks for your time.

Best,

Marco.

[Florian Augustin]

Hi Marco,

Thank you for reaching out to our technical support! This helps us a lot to find and fix issues in the product. My colleagues have looked at the error and will reach out to you shortly if they haven’t already.

To your other question regarding the steep transition in your model’s response when transitioning between different phases of the dynamics. In general, the model response will be non-differential at transition time if the change in model dynamics between the phases is discontinuous. This is also the case in my example above. However, the reactions are mass-action and the concentrations and changes in coefficients are small so that the curves still look relatively smooth visually.

Here is an idea how you could achieve smooth transitions using a sigmoid function. This approach might require some tweaking to get the transition period right, but maybe a parameter estimation (Fit program in the SimBiology Model Analyzer) could help if you have data.

Example:

First you would have to create a sigmoid-like function on your MATLAB path:

function transition = phase(time, transitionTime, timeScale)

% Smoothly transition between 0 and 1.

% * transitionTime specifies the time at which transition is 0.5.

% * timeScale specifies how rapid the transition between 0 and 1 is;

% the smaller the timeScale the faster the transition. To transition

% between 1 and 0, specify a negative value for timeScale.

transition = 1 ./ ( 1 + 10*exp(-(time-transitionTime)/timeScale));

end

Alternatively, you could also add a repeated assignment rule and a dummy parameter "transition" to your SimBiology model. Using the MATLAB function above, you can then implement smooth transitions in your SimBiology model as follows:

% Define model

model = sbiomodel("smoothed time-dependent reaction rates");

addspecies(model, "A", 10);

addparameter(model, "k1", 0.1);

addparameter(model, "k2", 1.0);

% Define reaction rate: here I am using rate1 = k1*A, rate2 = -k2 as an example.

reaction = addreaction(model, "A -> null");

reaction.ReactionRate = "phase(time, 3, -1)*k1*A - phase(time, 8, 1)*k2";

% [EDIT/UPDATE]: in order to get the plot shown below, you have to use

% reaction.ReactionRate = "phase(time, 3, -0.1)*k1*A - phase(time, 8, 0.1)*k2";

% Simulate model

configset = getconfigset(model);

configset.StopTime = 10;

configset.RuntimeOptions.StatesToLog = "A";

simData = sbiosimulate(model);

% Plot results

sbioplot(simData);

grid on

I hope this helps.

Best,

Florian

[mitpi_182]

Hello Florian,

Is there any chance you could show me how to do this using the Simbiology GUI? Specifically, how to integrate the sigmoid function in my model in order to use it to define the reaction rates. Many thanks for your time.

Best,

Marco.

[Florian Augustin]

Hi Marco,

attached is an sbproj file with the example above. The project is compatible with MATLAB version R2020b and later. I hope this helps you to get started.

Best,

Florian