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

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.)

Attachments:

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 :(

Attachments:

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?

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:

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, vy 
   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]"}]

Attachments:

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:

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]

Attachments:

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.