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]

Attachments:

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.

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?

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?

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

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?

Attachments:

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.

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

Attachments:

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.