How do I model stem cell growth using loops?

Questions
, ,
1

I am trying to model the growth of cancer stem cells using loops, but don’t have any experience modeling loops and this page wasn’t very helpful. The equations I have constructed thus far are as follows:

For[i = 1; c[0] = 1; c[1] = 1, i < 10, i++, 
     c[i_] := c[i - 1]*(1 + p1*v - (1 - p1) v - dr); Print[c[i]]]

For[i = 1; p[0] = 0, i < 10, i++, 
     p[i_] := c[i - 1] (1 - p1) v + p[i - 1] (1 + p2*v - (1 - p2) v - d); 
     Print[p[i]]]
    
For[i = 1; d[0] = 0, i < 10, i++, 
     d[i_] := p[i - 1] (1 - p2) v + t[i - 1] (1 - dr]); 
     Print[d[i]]]

I want to create a loop such that for i between 2 and n, I can determine the number of each type of cells (c=cancer stem cells, p=progenitor cells, d=terminally differentiated cell) there are symbolically. Once this is determined, I would also like to include a loop such that when the value of receptor occupancy, as defined as

receptor occupancy= ((c[i]*f1 + p[i]*f2 + d[i]*f3) r)/k

is less than than some value x1, then v is changed to vy1 and when it is less than some other value x2, the value of p1 is changed to p1y2.

Any help would be great!

2

Welcome to Wolfram Community Carolyn !

A few comments / questions:

  • You should carefully review the code you posted for typos. At least in the 3rd loop there is missing code in Subscript[dr, ] after coma. Note your code is formatted properly - you know it because you see syntax coloring. Click “Edit” in bottom right corner of your post to correct typos. But keep proper code formatting using the 1st button above or simply starting each code line with 4 spaces.

  • What do you mean “there are symbolically” ? That you do not specify values p1, dr, etc. numerically and want formulas in terms of these variables?

  • Is there any reference (article, Wikipedia, etc.) to these equations? Perhaps we could find a more efficient way to compute your results.

  • Keep in mind you can attach files to posts. Also this is very useful to read if you have not yet.

3

Hi Carolyn,

Below is an example of how a recursive system of definitions can be used for such a calculation. It makes use of the basic principle of the evaluation loop in Mathematica: When you evaluate an expression Mathematica continues to evaluate subsequent output until no built-in or user defined rule applies.

Using this, the intertwined recursive definitions provide a result for any index of any of the defined functions. In symbolic form it will of course use a lot of memory to get to high-index values. The earlier numerical values are substituted, the less memory required.

** Note that there is an issue in evaluation, because t[i] appears in the definition of d[i], but is not defined, it appears unevaluated in the results. **

In[1]:= (* Mathematica's evaluation engine handles recursion *)
(* Specific cases take precedence *)
(* Here is how we define and use c *)

In[2]:= c[0] = 1; c[1] = 1;

In[3]:= (* The = stores the value, so we don't have to calculate it \
again *)
c[i_] := c[i] = c[i - 1]*(1 + p1*v - (1 - p1) v - dr)

In[4]:= (* example *)
c[5]

Out[4]= (1 - dr - (1 - p1) v + p1 v)^4

In[5]:= (* and likewise for p *)

In[6]:= p[0] = 0;

In[7]:= p[i_] := 
 p[i] = c[i - 1] (1 - p1) v + p[i - 1] (1 + p2*v - (1 - p2) v - d)

In[8]:= (* example *)
p[5]

Out[8]= (1 - p1) v (1 - dr - (1 - p1) v + p1 v)^3 + (1 - 
    d - (1 - p2) v + 
    p2 v) ((1 - p1) v (1 - dr - (1 - p1) v + p1 v)^2 + (1 - 
       d - (1 - p2) v + 
       p2 v) ((1 - p1) v (1 - dr - (1 - p1) v + p1 v) + (1 - 
          d - (1 - p2) v + 
          p2 v) ((1 - p1) v + (1 - p1) v (1 - d - (1 - p2) v + p2 v))))

In[9]:= (* and for d *)

In[10]:= d[0] = 0;

In[11]:= (* I deleted the subscript typo *)
d[i_] := p[i - 1] (1 - p2) v + t[i - 1] (1 - dr)

In[12]:= (* example *)
d[11] // Simplify

Out[12]= (-1 + p1) (1 - 
    p2) v^2 (-(-1 + dr + v - 2 p1 v)^8 + (1 - 
       d + (-1 + 2 p2) v) ((-1 + dr + v - 
         2 p1 v)^7 + (1 - 
          d + (-1 + 2 p2) v) (-(-1 + dr + v - 2 p1 v)^6 + (1 - 
             d + (-1 + 2 p2) v) ((-1 + dr + v - 
               2 p1 v)^5 + (1 - 
                d + (-1 + 2 p2) v) (-(-1 + dr + v - 2 p1 v)^4 + (1 - 
                   d + (-1 + 2 p2) v) ((-1 + dr + v - 
                    2 p1 v)^3 + (1 - 
                    d + (-1 + 2 p2) v) (-(-1 + dr + v - 
                    2 p1 v)^2 + (-1 + d + v - 2 p2 v) (3 + d^2 - dr - 
                    4 v + 2 p1 v + 6 p2 v + v^2 - 4 p2 v^2 + 
                    4 p2^2 v^2 + d (-3 + (2 - 4 p2) v))))))))) - (-1 +
     dr) t[10]
4

Sam, thanks for your advice. Based on David’s advice on another post, I meant to remove all subscripts as they often cause problems, but must have missed that one. It is now corrected!

Yes, by “symbolically” that is what I mean. is there a more precise way of saying that?

There is no reference for these materials as it is a new proposed model for cell growth, one that I am hoping to investigate.

5

The t[i] is a type and should instead be d[i]. I changed this in the original question. Thanks for catching that. David, your code is much more eloquent than mine; thanks for the lesson recursive functions! I feel like I have a good handle on them but remain confused as to how I should incorporate the second rule, explained at the bottom of my post. Ideally at the end, I would like to be able to input different values for the variables and observe the behavior of the graphs as i changes.

6

These equations look like finite difference approximations to differential equations. That being the case, I think that writing them as such and using NDSolve or related functions is the way to do this. For example, ParametricNDSolve can generate a solutions set based on stated parameters, which can be investigated graphically. The last relationship you mention says that v and p1 are functions of the dependent variables. That can be included in the equations provided NDSolve. Discontinuity is OK.

I attach a fun example. It’s not parametric, but does include discontinuous algebraic functions of the dependent variables. A rocket is launched. The thrust occurs between tOn and Toff. At chuteTime a parachute is deployed, which significantly increases drag.

How about if we reformulate these as diffeq’s?

In[1]:= SetOptions[ParametricPlot, PlotTheme -> "Scientific"];

In[2]:= SetOptions[Plot, PlotTheme -> "Scientific"];

In[3]:= tImpact = 
  100; (* initialize - without this WhenEvent can't access it *)

In[4]:= (* rocket thrust *)
thrust[t_] := If[tOn <= t < tOff, f, 0]

In[5]:= (* drag coefficients without / with parachute *)
dragC[t_] := If[t < chuteTime, c1, c2]

In[6]:= (* values for simulation *)
values = {v0 -> 0.0001, theta -> 70 Degree, tOn -> 0, tOff -> 5, 
   f -> 150, chuteTime -> 40, c1 -> 0.0001, c2 -> 0.002, g -> 9.8};

In[7]:= equations = {
   {x''[t], y''[t]} == 
    thrust[t] Normalize[{x'[t], y'[t]}] - {0, g} - 
     dragC[t] Normalize[{x'[t], y'[t]}] Norm[{x'[t], y'[t]}]^2,
   {x'[0], y'[0]} == v0 {Cos[theta], Sin[theta]},
   {x[0], y[0]} == {0, 0},
   WhenEvent[y[t] < 0, {"StopIntegration", tImpact = t}]
   };

In[8]:= sol = 
  NDSolveValue[equations /. values, {x[t], y[t]}, {t, 0, 500}];

In[9]:= sol = 
  NDSolveValue[
   equations /. values, {x[t], y[t], x'[t], y'[t], x''[t], 
    y''[t]}, {t, 0, 500}];

In[10]:= tImpact

Out[10]= 90.9562

In[11]:= ParametricPlot[sol[[{1, 2}]], {t, 0, tImpact}, 
 AspectRatio -> Automatic, PlotRange -> All]

Rocket1_1.png
Rocket1_1.png360×133 2.09 KB

7

Ok so if I reformulate them as differential equations, then they can be solved using DSolve as follows:

DSolve[{c'[t] == c[t] (p1*v - (1 - p1) v - dr),
  p'[t] == c[t] (1 - p1) v + p[t] (p1v - (1 - p1) v - dr),
  d'[t] == p[t] (1 - p1) v + d[t] (1 - dr)}, {c, p, d}, t]

How do I proceed from here?

8

Here’s what I had in mind. It includes the discontinuous algebraics, although for my choice of parameter values they never switch.
I suggest you provide some “reasonable” parameter values, most of which are fixed, but a few for which you would like to see results for a range of values. (My values likely make no sense.)

In[1]:= (* function for occupancy *)
occupancy[c_, p_, d_] := r (c f1 + p f2 + d f3)/k

In[2]:= (* equations *)
eqs = {c'[t] == c[t] (p1[t]*v[t] - (1 - p1[t]) v[t] - dr), 
   p'[t] == 
    c[t] (1 - p1[t]) v[t] + p[t] (p1v - (1 - p1[t]) v[t] - dr), 
   d'[t] == p[t] (1 - p1[t]) v[t] + d[t] (1 - dr)};

In[3]:= (* initial conditions *)
initial = {c[0] == 1, p[0] == 0, d[0] == 0};

In[4]:= (* algebraic equations *)
(* values that vary with the independent parameter must be functions \
of it *)
algebraics = {v[t] == If[occupancy[c[t], p[t], d[t]] < x1, vy1, v0], 
   p1[t] == If[occupancy[c[t], p[t], d[t]] < x2, p1y2, p10]};

In[5]:= (* real numbers *)
values = {x1 -> .5, vy1 -> .5, v0 -> .3, x2 -> .6, p10 -> .3, 
   p1y2 -> .7, p10 -> .2, dr -> .5, p1v -> .4, f1 -> .5, f2 -> .3, 
   f3 -> .7, r -> .7, k -> 10};

In[6]:= (* check to make sure no undefined symbols *)
{eqs, initial, algebraics} /. values

Out[6]= {{Derivative[1][c][t] == 
   c[t] (-0.5 - (1 - p1[t]) v[t] + p1[t] v[t]), 
  Derivative[1][p][t] == 
   c[t] (1 - p1[t]) v[t] + p[t] (-0.1 - (1 - p1[t]) v[t]), 
  Derivative[1][d][t] == 0.5 d[t] + p[t] (1 - p1[t]) v[t]}, {c[0] == 
   1, p[0] == 0, 
  d[0] == 0}, {v[t] == 
   If[0.07 (0.5 c[t] + 0.7 d[t] + 0.3 p[t]) < 0.5, 0.5, 0.3], 
  p1[t] == If[0.07 (0.5 c[t] + 0.7 d[t] + 0.3 p[t]) < 0.6, 0.7, 0.3]}}

In[7]:= (* solve to f[t] set *)
{fc, fp, fd, fv, fp1} = NDSolveValue[
   {eqs, initial, algebraics} /. values,
   {c[t], p[t], d[t], v[t], p1[t]},
   {t, 0, 10}
   ];

In[8]:= Plot[{fc, fp, fd, fv, fp1}, {t, 0, 10}, 
 PlotLegends -> {"c[t]", "p[t]", "d[t]", "v[t]", "p1[t]"}]

8782cancermodel7.png
8782cancermodel7.png432×223 5.98 KB

9

Wow that’s amazing! That “If” will be really helpful for later projects. Can you explain what the command under “check to make sure no undefined symbols” is as far is the input, so that I know how to use it later? For example, if I didn’t have a code that had the equations, initial values, and algebraics outlined as you do, how would I check to make sure all my symbols are defined?

Also, I’ve seen “/.” used in other codes but I haven’t been able to find documentation as to what it does. Could you explain that as well?

10

.

11

Hi Carolyn, I think you saw incorrect code I posted, before I managed to correct it. Sorry for the confusion.

Regarding the check for undefined symbols – it is actually not necessary. But what I am checking is to see that the equations after the parameters are replaced with numerical values do not contain undefined values. They can contain the functions being solved for ( as some f[t] ) or other functions of independent or dependent variable which will resolve to numerical values. But if there are still symbolic parameters NDSolve will return the system unevaluated. It has to be able to get numerical values as it integrates. This can be seen in any case since these undefined symbols will be in the returned and partially evaluated result, but I often check first since it’s easier to read. I am just performing the Replace on the equation set and checking that it’s ready for NDSolve. (See below.)

/. is a short form infix operator for Replace.

In[13]:= Replace[x, x -> a]

Out[13]= a

Can be written as

In[15]:= x /. x -> a

Out[15]= a

Mathematica is a great symbolic algebra tool, and I prefer to work in symbols as long as convenient, and use a list of rules to substitute numerical values when I want numerical results, or for functions that require them, like NDSolve.

And since rules can be defined using patterns, and used within function definitions, they are often a convenient and readable way to define functions:

In[16]:= swap[z_] := z /. {x_, y_} -> {y, x}

In[17]:= swap[{a, b}]

Out[17]= {b, a}

It is worth noting that Mathematica is a rules engine. When you evaluate an expression, Mathematica’s evaluation engine looks to see if any rule may be applied. Here “rules” are those established by built-in functions or by user definitions made by assignments. If it does, it applies the rule, giving precedence to the most specific. (That’s at the heart of the recursive function definition.) After applying the rule, it looks again. And it keeps going until it obtains a result to which no rule applies. That is also why when you give a function and invalid argument it returns unevaluated. It looked, nothing applied, so the main evaluation loop terminates and gives you the first result to which no rule applied – Which in this case is your input.

You can easily change parameter values in the code above by changing the rule list values. But by using these methods, we can go further with NDSolve. We could define a function which has one or more of the parameter values as arguments, and uses Replace to put them in the equation set before using Replace for the parameter set as a whole. Then we have a function which can be passed an argument and generate a solution set. We could even use the solutions within our function definition to determine values at some time t, and return those values. Then we could, for example, plot some terminal concentrations as a function of one or more parameter values. Complicated function are often defined within a Module.

Mathematica is a bit of a steep learning curve, but once you have it, you have a Swiss Army Knife. It is really worth it. There is a resource I might have mentioned called The Virtual Book which can be found online, or by searching for it in the help system. It’s good reading.

Kind regards,
David

12

I thought it would be a good idea to use the command Manipulate so that I could observe the changes in the graph that correspond to changing the values of the constants. To do this, I need the symbolic solution to the system, which I found as follows:

equations = {c'[t] == c[t] (p1*v - (1 - p1) v - dr), 
  p'[t] == c[t] (1 - p1) v + p[t] (p1*v - (1 - p1) v - dr), 
  d'[t] == p[t] (1 - p1) v + d[t] (1 - dr)}


 DSolve[{equations, c[0] == 1, p[0] == 0, d[0] == 0}, {c[t], d[t], 
      p[t]}, {t}]
    ffc[t_] == E^(t (-dr + (-1 + 2 p1) v))
    ffp[t_] == 
     1/(-1 - v + 2 p1 v)^2 (-1 + p1)^2 v^2 (E^((1 - dr) t) - E^(
        t (-dr + (-1 + 2 p1) v)) - E^(t (-dr + (-1 + 2 p1) v)) t - 
        E^(t (-dr + (-1 + 2 p1) v)) t v + 
        2 E^(t (-dr + (-1 + 2 p1) v)) p1 t v)
    ffd[t_] == -E^(t (-dr + (-1 + 2 p1) v)) (-1 + p1) t v
  
  total[t_] == ffc[t] + ffp[t] + ffd[t]
   
 Manipulate[
     Plot[{ffc[t], ffp[t], ffd[t],total[t]}, {t, 0, 10}, 
      PlotLegends -> {"c[t]", "p[t]", "d[t]", "total cells"}], {{x1, 1, 
       "Occupany Level at which Division Speed Increases"}, 1, 
      10}, {{y1, 1, "Division Speed Increase Factor"}, 0, 
      10}, {{x2, 1, 
       "Occupancy Level at which Prop. Symmetric Division Increases"}, 1, 
      10}, {{y2, 1, "Division Proportion Increase Factor"}, 0, 
      10}, {{v0, 1, "Cell Divisions Per Unit Time"}, 0, 
      10}, {{p10, .5, 
       "Proportion of Cells Undergoing Symmetric Division"}, 0, 
      1}, {{dr, .2, "Cell Death Rate"}, 0, 
      1}, {{f1, 1, "CSC DLL4 Production"}, 0, 
      10}, {{f2, 1, "PC DLL4 Production"}, 0, 
      10}, {{f3, 1, "TDC DLL4 Production"}, 0, 
      10}, {{r, 1, "Receptors Per Cell"}, 0, 10}, {k, 0, 10}]

This is a very rough attempt and I know that I should include v[t] and p1[t] in the plot too, in addition to several other things, but I can’t seem to make the plot show any thing currently. Once I can get it to work I plan on tinkering with it. Anyexplanation as to why nothing shows up?

13

I thought it would be a good idea to use the command Manipulate so that I could observe the changes in the graph that correspond to changing the values of the constants. To do this, I need the symbolic solution to the system, which I found as follows:

equations = {c'[t] == c[t] (p1*v - (1 - p1) v - dr), 
  p'[t] == c[t] (1 - p1) v + p[t] (p1*v - (1 - p1) v - dr), 
  d'[t] == p[t] (1 - p1) v + d[t] (1 - dr)}


 DSolve[{equations, c[0] == 1, p[0] == 0, d[0] == 0}, {c[t], d[t], 
      p[t]}, {t}]
    ffc[t_] == E^(t (-dr + (-1 + 2 p1) v))
    ffp[t_] == 
     1/(-1 - v + 2 p1 v)^2 (-1 + p1)^2 v^2 (E^((1 - dr) t) - E^(
        t (-dr + (-1 + 2 p1) v)) - E^(t (-dr + (-1 + 2 p1) v)) t - 
        E^(t (-dr + (-1 + 2 p1) v)) t v + 
        2 E^(t (-dr + (-1 + 2 p1) v)) p1 t v)
    ffd[t_] == -E^(t (-dr + (-1 + 2 p1) v)) (-1 + p1) t v
  
  total[t_] == ffc[t] + ffp[t] + ffd[t]
   
 Manipulate[
     Plot[{ffc[t], ffp[t], ffd[t],total[t]}, {t, 0, 10}, 
      PlotLegends -> {"c[t]", "p[t]", "d[t]", "total cells"}], {{x1, 1, 
       "Occupany Level at which Division Speed Increases"}, 1, 
      10}, {{y1, 1, "Division Speed Increase Factor"}, 0, 
      10}, {{x2, 1, 
       "Occupancy Level at which Prop. Symmetric Division Increases"}, 1, 
      10}, {{y2, 1, "Division Proportion Increase Factor"}, 0, 
      10}, {{v0, 1, "Cell Divisions Per Unit Time"}, 0, 
      10}, {{p10, .5, 
       "Proportion of Cells Undergoing Symmetric Division"}, 0, 
      1}, {{dr, .2, "Cell Death Rate"}, 0, 
      1}, {{f1, 1, "CSC DLL4 Production"}, 0, 
      10}, {{f2, 1, "PC DLL4 Production"}, 0, 
      10}, {{f3, 1, "TDC DLL4 Production"}, 0, 
      10}, {{r, 1, "Receptors Per Cell"}, 0, 10}, {k, 0, 10}]

This is a very rough attempt and I know that I should include v[t] and p1[t] in the plot too, in addition to several other things, but I can’t seem to make the plot show any thing currently. Once I can get it to work I plan on tinkering with it. Anyexplanation as to why nothing shows up?

14

In the lines after DSolve, you are using == rather than = to define functions. For example:

ffc[t_] == E^(t (-dr + (-1 + 2 p1) v))

Should be

ffc[t_] = E^(t (-dr + (-1 + 2 p1) v)) ( = is short for Set )

Or most often

ffc[t_] := E^(t (-dr + (-1 + 2 p1) v)) ( := is short for SetDelayed)

The difference is that = will evaluate the right hand side before using it as a template to define the function. := defers the evaluation of the RHS until the function is invoked. Unless you know you need to evaluate the RHS, := (SetDelayed) is the usual choice.

== asserts equality in an equation, or tests for equality in a conditional. So none of these actually defined a function.

When you get blank plots there are two good steps to take. The first, and most useful, is to take the expression being plotted, by itself, and use a rule replace to evaluate it at a numerical value. For example, if the independent variable in the expression is x, evaluate expression/.x->10, or some real number in the domain. You need it to return a numerical result, because that’s what Plot needs. If it doesn’t, you will see what it is in the expression that is preventing evaluation from yielding a number.The second is to make a simple-as-possible plot for the expression, meaning without Manipulate.

So, first rewrite your functions with :=

Then try evaluating the expression you are plotting with a number.

Best,

David

15

I should add that an expression intended for manipulate will have one or more variables undefined, so manipulate can make the substitutions. These also need numerical values during Plotting, but they will be provided by Manipulate. So for such an expression the Replace test indicates success if you evaluate the expression with a numerical independent variable and the only undefined symbols remaining are those which will be parameters in Manipulate.

16

I performed the two tests you advised, and was able to plot each of the individual equations without Manipulate. When I applied Manipulate to one of them (I chose p(t)) and only established one constant to be manipulated, the graph sadly remained blank :frowning:

17

Hi Carolyn,

There was an additional problem, and it’s somewhat advanced. It has to do with how Plot and Manipulate localize variables. The short solution is that all of the parameters being manipulated should appear explicitly within the Plot expression. I attach this as a notebook.

You should also study the way the result of DSolve is being used to define the functions we want from its output. DSolve outputs rules, and we use these to generate expressions, which are themselves used as the templates for function definitions. And so these templates allow us to define the functions as having their parameters as arguments, as well as the independent variable t.

I’m not usually this helpful ;-). But it is 95 degrees heading for 100 here in the state of Washington. We prefer 65-75. So I’m indoors amusing myself with Mathematica.

Best, David

18

And here is one more version which I think is better. Rather than build functions with all the parameters exposed as arguments, we build functions of t, which still have undefined variables internally. We expose these to manipulate not as explicit arguments, but rather we use dummy parameters in manipulate, so manipulate can see them, but then we use a Replace with a Rule in the expression being plotted. This replace causes the “hidden” parameters in the expression to be replaced by the numerical values generated for the dummy parameters by manipulate. I think this is better because the functions are not so cluttered up with a long argument list. They can always be used with specific values for the undefined variable by using a similar Replace in other applications. (Note that I did remove a lot of unused parameters from manipulate.)

19

When using Manipuate, I simplified the equation by removing the occupancy function from the mix and allowing the v[t] and p1[t] be constants rather than functions. Is there a way to reintroduce these easily? I originally thought I could use the code you constructed and then a replace command that would replace v->v[t] and p1->p1[t] but I can’t seem to find something that will substitute just portions of an equation. Even after this replacement, I am concerned the code wouldn’t work in Manipulate since based off what you have shown me, it seems to be very sensitive to the type of problem I am working on.

EDIT:
Using /. (thanks for teaching me about that) I was able to replace what I needed so that the equations read how I want them to (i.e. include the functions p1[t] and v[t] as you define in your earlier code. The issue with using Manipulate to graph them persists however.

20

One way to handle occupancy[] is to include it in the equation set as occupancy[t] == expression, just like the other equations, rather than as a defined function. I had been thinking this earlier because it would permit occupancy to be included in the the list of functions or rules which DSolve or NDSolve returns. Then it could be graphed with the rest so its behavior can be observed.

Regarding the Manipulate problem. I think it’s reasonable to say that it is sensitive not to the type of problem, but to exact syntax. I would need to see what you have done to be able to comment. (If you prefer not to post it in this forum, you can also contact me through email – look up my profile here, and follow to my web site.)

Best,
D