How can I write this system of differential equations so MMA can solve it?

Questions
,
1

I am trying to solve a system of nonlinear differential equations, which I set up with the following code:

B = {{(Subscript[p, 0] - Subscript[q, 0])*(Subscript[v, 
       0]/(1 + Subscript[\[Beta], 0] ((Subscript[x, 2][t])^2))) - 
    Subscript[d, 0], 0, 
   0}, {(1 - Subscript[p, 0] + Subscript[q, 0])*(Subscript[v, 
      0]/(1 + Subscript[\[Beta], 
         0] ((Subscript[x, 2][t])^2))), (Subscript[p, 1] - Subscript[
       q, 1])*(Subscript[v, 
       1]/(1 + Subscript[\[Beta], 1] ((Subscript[x, 2][t])^2))) - 
    Subscript[d, 1], 
   0}, {0, (1 - Subscript[p, 1] + Subscript[q, 1])*(Subscript[v, 
      1]/(1 + Subscript[\[Beta], 
         1] ((Subscript[x, 2][t])^2))), -Subscript[d, 2]}}
    

    X[t_] = {Subscript[x, 0][t], Subscript[x, 1][t], Subscript[x, 2][t]};
    systemB = X'[t] == B.X[t]

solB = DSolve[
  systemB, {Subscript[x, 0], Subscript[x, 1], Subscript[x, 2]}, t]

This code returns the exact input, which from my research means that it cannot be solved or takes too long to solve. I spoke to someone who wrote a paper including this equation and he said Mathematica should be able to solve it. Am I writing it incorrectly? Even with the initial conditions

Subscript[x, 1][ 0] == 0, Subscript[x, 2][ 0] == 0

I can’t seem to get it to work.

Any help would be great!

2

Hi Carolyn,

Two suggestions: 1) Do not use subscripted variables. 2) Please see How to post for information on how to put your code in a code block. That way we can copy-paste it into Mathematica and try it out. (When I try that with your definition for B I get a syntax error.)

Best Regards,
David

3

The code block is good. But the subscripts are still likely an issue. The problem is that Mathematica wants the Head for variables to be Symbol, but the Head for a subscripted character is Subscript. There are ways to make this work, but most of us avoid them like the plague.

In[9]:= Head[p0]

Out[9]= Symbol

In[10]:= Head[Subscript[p, 0]]

Out[10]= Subscript
4

I removed all subscripts so that my equation reads

DSolve[{x'[t] == ((a - b)*(c/(1 + n ((z[t])^2))) - j)*x[t], 
  y'[t] == ((1 - a + b)*(c/(1 + n ((z[t])^2))))*
     x[t] + ((f - g)*(h/(1 + r ((z[t])^2))) - k)*y[t], 
  z[t] == ((1 - f + g)*(h/(1 + r ((z[t])^2))))*y[t] + (-m)*z[t], 
  y[0] == 0, z[ 0] == 0, x[ 0] == 1}, {x, y, z}, t]

but all that is returned is the input. Any additional advice?

5

I did the same – see below. And I got the same result. Either we are both making an error, or Mathematica does not find a symbolic solution to this system. NDSolve might provide a numerical solution for specific parameter values and initial conditions. Or perhaps someone familiar with this system could comment. :frowning:

In[1]:= B = {{(p0 - q0)*(v0/(1 + \[Beta]0 ((x2[t])^2))) - d0, 0, 
   0}, {(1 - p0 + 
      q0)*(v0/(1 + \[Beta]0 ((x2[t])^2))), (p1 - 
       q1)*(v1/(1 + \[Beta]1 ((x2[t])^2))) - d1, 
   0}, {0, (1 - p1 + q1)*(v1/(1 + \[Beta]1 ((x2[t])^2))), -d2}}


Out[1]= {{-d0 + ((p0 - q0) v0)/(1 + \[Beta]0 x2[t]^2), 0, 
  0}, {((1 - p0 + q0) v0)/(
  1 + \[Beta]0 x2[t]^2), -d1 + ((p1 - q1) v1)/(1 + \[Beta]1 x2[t]^2), 
  0}, {0, ((1 - p1 + q1) v1)/(1 + \[Beta]1 x2[t]^2), -d2}}

In[2]:= X[t_] = {x0[t], x1[t], x2[t]}

Out[2]= {x0[t], x1[t], x2[t]}

In[3]:= systemB = X'[t] == B.X[t]

Out[3]= {Derivative[1][x0][t], Derivative[1][x1][t], 
  Derivative[1][x2][
   t]} == {x0[
    t] (-d0 + ((p0 - q0) v0)/(
     1 + \[Beta]0 x2[t]^2)), ((1 - p0 + q0) v0 x0[t])/(
   1 + \[Beta]0 x2[t]^2) + 
   x1[t] (-d1 + ((p1 - q1) v1)/(1 + \[Beta]1 x2[t]^2)), -d2 x2[
     t] + ((1 - p1 + q1) v1 x1[t])/(1 + \[Beta]1 x2[t]^2)}

In[4]:= solB = DSolve[systemB, {x0, x1, x2}, t]

Out[4]= DSolve[{Derivative[1][x0][t], Derivative[1][x1][t], 
   Derivative[1][x2][
    t]} == {x0[
     t] (-d0 + ((p0 - q0) v0)/(
      1 + \[Beta]0 x2[t]^2)), ((1 - p0 + q0) v0 x0[t])/(
    1 + \[Beta]0 x2[t]^2) + 
    x1[t] (-d1 + ((p1 - q1) v1)/(1 + \[Beta]1 x2[t]^2)), -d2 x2[
      t] + ((1 - p1 + q1) v1 x1[t])/(1 + \[Beta]1 x2[t]^2)}, {x0, x1, 
  x2}, t]
6

They system is from this paper on cancer stem cells, in case that helps anyone with solving the equation (equation S4 in the supplemental materials of the paper).

Since Mathematica doesn’t seem able to solve this, do you know another program that might be able to?

7

I entered equations directly from the paper, and DSolve finds a solution. However, the symbolic solution violates the initial condition for x2[t]. I suspect some error on my part, but have not found it. I attach the notebook in the hope that other eyes will see what mine have not.
Best,
David

8

And on futher study, it appears this system of equations does not have an initial solution when initial conditions are defined for the three dependent variables:

In[42]:= initialStep = (eqs /. {x0[t] -> x00, x1[t] -> x10, 
     x2[t] -> x20})[[1 ;; 3]]

Out[42]= {Derivative[1][x0][t] == -d0 x00 + (p0 - q0) v0 x00, 
 Derivative[1][x1][t] == (1 - p0 + q0) v0 x00 - 
   d1 x10 + (p1 - q1) v1 x10, 
 Derivative[1][x2][t] == (1 - p1 + q1) v1 x10 - d2 x20}

In[43]:= initialStepEqs = 
 initialStep /. {x0'[t] -> dxo, x1'[t] -> dx1, x2'[t] -> dx2}

Out[43]= {dxo == -d0 x00 + (p0 - q0) v0 x00, 
 dx1 == (1 - p0 + q0) v0 x00 - d1 x10 + (p1 - q1) v1 x10, 
 dx2 == (1 - p1 + q1) v1 x10 - d2 x20}

In[44]:= Solve[initialStepEqs, {dx0, dx1, dx2}]

Out[44]= {}
9

The equations you used are the more simple versions without feedback loops. I solved those using DSolve as well, however the more complicated one, where v0, v1, p0, p1, q0, and q1 are replaced by nonlinear functions (shown by equations 2 and 3) was where I ran into trouble.

I really appreciate all you help!

10

Yes – and I see that in the simple version the the equation for x0[t] stands alone, and can be used then in the solution for the next two. I did not look closely enough to see the difference in your second set. I suspect a good course is to use NDSolve to get a numerical solution. Do you have numerical values for the parameters and initial conditions?

11

Ideally I was hoping to use the symbolic solution, but I can find numerical values to work on it. You’ve been incredibly helpful. If you have more time I have another question posted here about loops. I’m new to mathematica so the learning curve is very steep.

12

Here is an example of NDSolveValue being used to obtain a numerical solution to the simpler equations found in the paper:

In[9]:= (* equations and initial conditions *)
eqs = {
   x0'[t] == (p0 - q0) v0 x0[t] - d0 x0[t],
   x1'[t] == (1 - p0 + q0) v0 x0[t] + (p1 - q1) v1 x1[t] - d1 x1[t],
   x2'[t] == (1 - p1 + q1) v1 x1[t] - d2 x2[t],
   x0[0] == x00, x1[0] == x10, x2[0] == x20
   };

(* made up parameter values *)
par = {q0 -> .5, q1 -> .5, v0 -> 1, v1 -> 1, x00 -> 100, x10 -> 0, 
   x20 -> 0, p0 -> .5, p1 -> .4, p2 -> .5, d0 -> .2, d1 -> .3, 
   d2 -> .4};

In[11]:= eqs /. par

Out[11]= {Derivative[1][x0][t] == 0. - 0.2 x0[t], 
 Derivative[1][x1][t] == 1. x0[t] - 0.4 x1[t], 
 Derivative[1][x2][t] == 1.1 x1[t] - 0.4 x2[t], x0[0] == 100, 
 x1[0] == 0, x2[0] == 0}

In[15]:= {fx0, fx1, fx2} = 
  NDSolveValue[eqs /. par, {x0, x1, x2}, {t, 0, 100}];

In[14]:= (* plot x0, x1, x2 with numerical values *)
LogPlot[{fx0[t], fx1[t], fx2[t]}, {t, 0, 100}, PlotRange -> All, 
 PlotLegends -> {"x0", "x1", "x2"}]

cancermodel5.png
cancermodel5.png360×213 3.35 KB

13

Hi Carolyn, I saw your loops. At first glance you appeared to be doing what NDSolve can do. (And your back to subscripts! Take my word for it – they are Trouble with a capital T.) ;-}

I’ll look at your more complicated set with NDSolve.

D

14

And here is a numerical solution for your nonlinear system, with made up parameter values. When forming systemB, I used Thread so equations would be formatted individually, so they could be joined with the initial conditions. I also attach the notebook.

In[1]:= B = {{(p0 - q0)*(v0/(1 + \[Beta]0 ((x2[t])^2))) - d0, 0, 
   0}, {(1 - p0 + 
      q0)*(v0/(1 + \[Beta]0 ((x2[t])^2))), (p1 - 
       q1)*(v1/(1 + \[Beta]1 ((x2[t])^2))) - d1, 
   0}, {0, (1 - p1 + q1)*(v1/(1 + \[Beta]1 ((x2[t])^2))), -d2}}

Out[1]= {{-d0 + ((p0 - q0) v0)/(1 + \[Beta]0 x2[t]^2), 0, 
  0}, {((1 - p0 + q0) v0)/(
  1 + \[Beta]0 x2[t]^2), -d1 + ((p1 - q1) v1)/(1 + \[Beta]1 x2[t]^2), 
  0}, {0, ((1 - p1 + q1) v1)/(1 + \[Beta]1 x2[t]^2), -d2}}

In[2]:= initB = {x0[0] == x00, x1[0] == x10, x2[0] == x20};

In[12]:= X[t_] = {x0[t], x1[t], x2[t]};

In[13]:= systemB = Thread[X'[t] == B.X[t]];

In[14]:= equations = Join[systemB, initB];

In[7]:= par = {q0 -> .5, q1 -> .5, v0 -> 1, v1 -> 1, x00 -> 100, 
   x10 -> 0, x20 -> 0, p0 -> .5, p1 -> .4, p2 -> .5, d0 -> .2, 
   d1 -> .3, d2 -> .4, \[Beta]0 -> 1, \[Beta]1 -> 1};

In[8]:= equations /. par

Out[8]= {Derivative[1][x0][t] == -0.2 x0[t], 
 Derivative[1][x1][t] == (1. x0[t])/(1 + x2[t]^2) + 
   x1[t] (-0.3 - 0.1/(1 + x2[t]^2)), 
 Derivative[1][x2][t] == -0.4 x2[t] + (1.1 x1[t])/(1 + x2[t]^2), 
 x0[0] == 100, x1[0] == 0, x2[0] == 0}

In[15]:= {fx0, fx1, fx2} = 
  NDSolveValue[equations /. par, {x0, x1, x2}, {t, 0, 100}];

In[17]:= LogPlot[{fx0[t], fx1[t], fx2[t]}, {t, 0, 100}, 
 PlotLegends -> {"x0", "x1", "x2"}, PlotTheme -> "Scientific", 
 FrameLabel -> {"Time(s)", "Concentration"}]

cancermodel6.png
cancermodel6.png417×229 5.62 KB

15

Thank you!

16

Hello Carolyn and David,

This looks like a really tough problem. I just saw this today and took a look at the paper, so I have not tried to solve it yet.

I have to say that if Mathematica cannot solve a differential equation symbolically, then you would be hard pressed to find another program that could. Not trying to say that Mathematica is the best thing ever (although I do think it is quite impressive), it is just that not many programs are good at symbolic manipulations. I recall sitting through a presentation about a MAPLE tool-bench that couples to MATLAB. They were discussing symbolic solution of differential equations in the context of space probe design or the like. You may want to check that. Unfortunately I do not have references handy.

I think part of the reason your problem becomes so hard is that when you introduce the functions evaluated at a different time in the right-hand side (x2(t-tau)), then strictly speaking you probably no longer have a differential equation (not any equation with derivatives in it is a differential equation). I will try to take a look when I have a chance, but it looks like a symbolic solution would be very hard.

I agree with David on the subscripts. They are nice to display results in elegant forms. I have spent hours chasing illogical results which at the end turn out to be a mess-up with a subscripted variable.

A separate question for David, did you use to live in Pittsburgh? It just so happens that I know a David Keith…

Best,

OL.

17

Hi Otto,

Carolyn and I have continued this discussion off-line. And her model is in fact even more complicated than that in the paper she referenced. Mathematica can of course solve systems of linear ODEs. And it does provide solutions to quite a few nonlinear ODEs. But in my experience, most nonlinear ODE’s do not have closed form solutions, Numerical methods are required. I am not surprised to find it does not solve a nonlinear system. (That is not to say there is no closed form solution, but I am lazy enough to accept its null response for evidence of that.) As to Pittsburgh – It must be another me – I’ve been a West Coast guy all my life, currently in Vancouver WA.

Kind regards,
David

18

Dave,

If it is more complicated that in the paper, then it is really complicated. I would settle for a good numerical solution.

I am not sure you understood my point: I was questioning whether the thing was a differential equation at all (linear or not). Once the behavior of your system depends on the value of some quantity some finite time ago, you lose locality in time. It may not be a differential equation even if it looks like one. For example, dx/dt = x(x(t)) is not a differential equation.

That may seem like a pedantic formality, but the reason I point it out is that if the thing is really not a differential equation, then most likely the approaches built into Mathematica do not work. Of course, there is always hope for a numerical approach.

As for Pittsburgh, what can I say, know I know two David Keiths. I am sure Vancouver is no worse…

OL

19

The supplementary information from the paper, points to delay differential equations, NDSolve has this capability, although you must specify the delay explicitly (See the help section). The delays are a convenient way to capture processes that take place outside of the scope of the normal evolution process. In theory you could replace these with a full model of the delayed process, however you would have to solve your model over multiple time-scales.

20

Hi Otto and Pat, and thanks to both of you.

Otto; Sorry I didn’t comment on (t-tau) – you are of course quite correct – it takes the equation into troubled waters. I skipped over it no doubt because the form CBro was using did not include the delay (another papers did not) – but it was still a nonlinear system for which I suspected no closed form solution exists. I have to admit that I did not read the material in depth, thinking to comment only on the Mathematica aspects.

And thanks, Pat – I have not used delays in diffeqs before. I knew generally that the capability exists, but I will look into it.

I don’t know if Carolyn is still following this discussion, but I hope she rejoins it.

Best, David