Time-dependent Growth Rate of Species in Model

Question

mitpi_182 2022-6-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1741355-time-dependent-growth-rate-of-species-in-model

编辑： Florian Augustin 2022-7-29

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Florian Augustin 2022-6-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1741355-time-dependent-growth-rate-of-species-in-model#answer_986530

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)

Plot showing simulation results for example.

Best,

Florian

9 个评论
显示 7更早的评论隐藏 7更早的评论

Florian Augustin 2022-6-25

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:

Simulation Settings in SimBiology Model Builder

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:

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.
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 2022-7-25

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 2022-7-27

编辑：Florian Augustin 2022-7-29

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 2022-7-28

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 2022-7-29

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

请先登录，再进行评论。