Plotting a phase portrait of a reduced system with 4 differential equations

Question

I have a system of differential equations that I have reduced to this:

\begin{equation} \small \frac{db}{dt}= \beta_1 (1-b(t)-c(t)) b(t) + (1-\beta_1) \beta_3 (1-b(t)-c(t)) d(t) - \phi_{w1} b(t) c(t) - \phi_{b1} b(t) e(t) - \mu_B b(t) - \gamma_1 b(t) \end{equation} \begin{equation} \small \frac{dc}{dt}= \beta_2 (1-b(t)-c(t)) c(t) + (1-\beta_2) \beta_4 (1-b(t)-c(t)) e(t) + \phi_{w1} b(t) c(t) + \phi_{b1} b(t) e(t) - \mu_C c(t) - \gamma_2 c(t) \end{equation} \begin{equation} \small \frac{dd}{dt}= \beta_3 (1-d(t)-e(t)) d(t) + (1-\beta_3) \beta_1 (1-d(t)-e(t)) b(t) - \phi_{w2} d(t) e(t) - \phi_{b2} d(t) c(t) - \mu_D d(t) - \gamma_3 d(t) \end{equation} \begin{equation} \small \frac{de}{dt}= \beta_4 (1-d(t)-e(t)) e(t) + (1-\beta_4) \beta_2 (1-d(t)-e(t)) c(t) + \phi_{w2} d(t) e(t) + \phi_{b2} d(t) c(t) - \mu_E e(t) - \gamma_4 e(t) \end{equation}

Here the code for the equations as input in mathematica:

    ode1 = b'[t] == 
   Subscript[\[Beta], 1]*(1 \[Minus] b[t] \[Minus] c[t])*
     b[t] + (1 \[Minus] Subscript[\[Beta], 1])*
     Subscript[\[Beta], 3]*(1 \[Minus] b[t] \[Minus] c[t])*
     d[t] \[Minus] Subscript[\[Phi], w1]*b[t]*c[t] \[Minus] 
    Subscript[\[Phi], b1]*b[t]*e[t] \[Minus] 
    Subscript[\[Mu], B]*b[t] \[Minus] Subscript[\[Gamma], 1]*b[t];
ode2 = c'[t] == 
   Subscript[\[Beta], 2]*(1 \[Minus] b[t] \[Minus] c[t])*
     c[t] + (1 \[Minus] Subscript[\[Beta], 2])*
     Subscript[\[Beta], 4]*(1 \[Minus] b[t] \[Minus] c[t])*e[t] + 
    Subscript[\[Phi], w1]*b[t]*c[t] + 
    Subscript[\[Phi], b1]*b[t]*e[t] \[Minus] 
    Subscript[\[Mu], C]*c[t] \[Minus] Subscript[\[Gamma], 2]*c[t];
ode3 = d'[t] == 
   Subscript[\[Beta], 3] *(1 \[Minus] d[t] \[Minus] e[t])*
     d[t] + (1 \[Minus] Subscript[\[Beta], 3])*
     Subscript[\[Beta], 1]*(1 \[Minus] d[t] \[Minus] e[t])*
     b[t] \[Minus] Subscript[\[Phi], w2]*d[t]*e[t] \[Minus] 
    Subscript[\[Phi], b2]*d[t]*c[t] \[Minus] 
    Subscript[\[Mu], D]*d[t] \[Minus] Subscript[\[Gamma], 3]*d[t];
ode4 = e'[t] == 
   Subscript[\[Beta], 4] *(1 \[Minus] d[t] \[Minus] e[t])*
     e[t] + (1 \[Minus] Subscript[\[Beta], 4])*
     Subscript[\[Beta], 2]*(1 \[Minus] d[t] \[Minus] e[t])*c[t] + 
    Subscript[\[Phi], w2]*d[t]*e[t] + 
    Subscript[\[Phi], b2]*d[t]*c[t] \[Minus] 
    Subscript[\[Mu], E]*e[t] \[Minus] Subscript[\[Gamma], 4]*e[t];

Would the code to solve the system in mathematica be something like this?

  b[0] = 100
  c[0] = 500
  d[0] = 100
  e[0] = 400
  (fb, fc, fd, fe} = 
 NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == k1, c[0] == k2, 
    d[0] == k3, e[0] == k4}, {b, c, d, e}, {t, 0, 1000}].

I was trying to generate a visualisation of the Phase portrait of this system in mathematica.

Here b, c, d and e are positive real variables (population). Variable t goes from 0 to 1000. All the rest of parameters are all decimal numbers taking values from 0 to 1. We can assume constant prameters $\mu$ = 0.06 and constant parameters $\gamma$=0.01 in all the cases.

How can I a make a good visualisation of phase portraits of this system. For example, a 3D plot.

Can anyone help me?

Answer 1

It's a good try pyring. But in order to solve first-order ODEs numerically, you need two things: (1) Supply numeric values to all constants, (2) Supply numeric values for initial conditions. You have 16 constants and four equations so that's 20 values you need to assign. In the code below I chose random values between 0 and 1. The important thing is to just get it running without syntax errors or run-time errors. Later you can adjust it to your needs. In the code below I plotted individually the four functions written as fb[t], fc[t],fd[t] and fe[t]. Then I chose three of them fb[t], fc[t] and fd[t] to construct a 3D phase portrait which appears as a simple line. Maybe with other choices of initial conditions and constants you can achieve something more interesting. Just cut this code and paste it into your notebook and it should run without errors. Then you can experiment with changing the initial conditions and constants.

 (* 
     set up constants and ODEs
    *)
    Subscript[\[Mu], B] = 0.5;
     Subscript[\[Beta], 1] = 0.1;
    Subscript[\[Beta], 3] = 0.2;
    Subscript[\[Gamma], 1] = 0.7;
    Subscript[\[Gamma], 2] = 0.221;
    Subscript[\[Phi], b1] = 0.9;
    Subscript[\[Phi], w1] = 0.05;
    Subscript[\[Phi], w2] = 0.08;
    Subscript[\[Beta], 2] = 0.8;
    Subscript[\[Beta], 4] = 0.2;
    Subscript[\[Mu], C] = 0.5;
     Subscript[\[Mu], E] = 0.01;
    Subscript[\[Phi], b2] = 0.6;
    Subscript[\[Gamma], 3] = 0.109;
    Subscript[\[Mu], D] = .07;
    Subscript[\[Gamma], 4] = 0.219;
    ode1 = Derivative[1][b][t] == -(b[t] Subscript[\[Mu], B]) + 
        b[t] Subscript[\[Beta], 
         1] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta], 
           1]) Subscript[\[Beta], 3] d[t] (-b[t] - c[t] + 1) - 
        b[t] Subscript[\[Gamma], 1] - b[t] e[t] Subscript[\[Phi], b1] - 
        b[t] c[t] Subscript[\[Phi], w1];
    ode2 = Derivative[1][c][t] == 
       Subscript[\[Beta], 2]c[t] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta],2]) Subscript[\[Beta], 4] e[t] (-b[t] - c[t] + 1) + 
        b[t] e[t] Subscript[\[Phi], b1] + 
        b[t] c[t] Subscript[\[Phi], w1] - Subscript[\[Gamma], 2] c[t] - 
        c[t] Subscript[\[Mu], C];
    ode3 = Derivative[1][d][t] == 
       b[t] Subscript[\[Beta], 
         1] (1 - Subscript[\[Beta], 3]) (-d[t] - e[t] + 1) - 
        d[t] c[t] Subscript[\[Phi], b2] - Subscript[\[Gamma], 3] d[t] - 
        d[t] Subscript[\[Mu], D] + 
        Subscript[\[Beta], 3] d[t] (-d[t] - e[t] + 1) - 
        d[t] e[t] Subscript[\[Phi], w2];
    ode4 = Derivative[1][e][t] == 
       d[t] c[t] Subscript[\[Phi], b2] + 
        Subscript[\[Beta], 
         2] (1 - Subscript[\[Beta], 4]) c[t] (-d[t] - e[t] + 1) + 
        Subscript[\[Beta], 4] e[t] (-d[t] - e[t] + 1) + 
        d[t] e[t] Subscript[\[Phi], w2] - Subscript[\[Gamma], 4] e[t] - 
        e[t] Subscript[\[Mu], E];
    (*
     run system with some arbitrary initial conditions
    *)
    {fb, fc, fd, fe} = 
      NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == 0.5, c[0] == -0.2, 
        d[0] == 1.1, e[0] == 0.5}, {b, c, d, e}, {t, 0, 1000}];
    (*
     plot each function separately
    *)
    
    Plot[{fb[t], fc[t], fd[t], fe[t]}, {t, 0, 1000}, 
     PlotLabel -> Style["fb[t],fc[t],fd[t],fe[t]", 16]]
    (*
     choose {fb[t],fc[t],fd[t]} and plot a 3D phase portrait of these \
    fuctions
    *)
    ParametricPlot3D[{fb[t], fc[t], fd[t]}, {t, 0, 100}, 
     BoxRatios -> {1, 1, 1}, 
     PlotLabel -> Style["Phase portrait of {fb[t],fc[t],fd[t]}", 16]]