3
$\begingroup$

Welcome to all, I try to apply most suggestion codes of a plot bifurcation diagram but I could not, I got an error. Is there anyone who can help me to understand how I can plot a bifurcation diagram for the following system:

r = 0.431201; β1 = 2.99*10^-6; β2 = 1;
α1 = 0.44257; α2 = 0.4 ; α3 = 2.99*10^-6;
k1 = 0.11; k2 = 0.99; c1 = 0.33; c2 = 0.66;
sys3 = {
  N1'[t] == r N1[t] (1 - β1 N1[t]) - β2 N1[t] T[t] + c1 V[t] N1[t],
  T'[t] == α1 T[t] (1 - α2 T[t]) + α3 N1[t] T[t] - c2 V[t] T[t],
  V'[t] == k1 - k2 V[t]};
con3 = {N1[0] == 1, T[0] == 1, V[0] == 2};
s3 = NDSolve[{sys3, con3}, {N1, T, V}, {t, 0, 360}]

here, I want to plot the bifurction diagram of the system, also, the @zhk solve my problem for parametric solution code, but stil my qustion does not had a complete solution.

$\endgroup$

marked as duplicate by zhk, MarcoB, bbgodfrey differential-equations Mar 12 at 1:41

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • $\begingroup$ What do you want on the x and y axes of this bifurcation diagram? Also, could you say what the model represents? $\endgroup$ – Chris K Mar 4 at 7:07
  • $\begingroup$ I want to see the effect of equation 3 on the behaviour of N1and T this model describes the behaviour of normal and tumor cells. $\endgroup$ – sana alharbi Mar 4 at 9:41
  • 1
    $\begingroup$ Thanks, but I'm not sure what you mean by "equation 3". Is there one particular parameter of interest? $\endgroup$ – Chris K Mar 4 at 9:45
  • $\begingroup$ can study more than one parameter, actually, the behaviour of the model has affected by \beta_2, c_1, c_2 , \alhpha _3 and k2 but I do not if I can examine all these parameters where already I know that the values of them $\endgroup$ – sana alharbi Mar 4 at 9:51
  • $\begingroup$ @sanaalharbi How about this? ParametricPlot3D[Evaluate@{N1[t], T[t], V[t]} /. s3, {t, 0, 360}, PlotRange -> {{0, 3}, {0, 3}, {0, 3}}] $\endgroup$ – zhk Mar 4 at 16:10
3
$\begingroup$

I think you are looking for something like this

r = 0.431201; \[Beta]1 = 2.99*10^-6;\[Alpha]1 = 0.44257; \[Alpha]2 = 0.4;

Manipulate[

eq1 = x'[t] == r x[t] (1 - \[Beta]1 x[t]) - \[Beta]2 x[t] y[t] + c1 z[t] x[t];
eq2 = y'[t] == \[Alpha]1 y[t] (1 - \[Alpha]2 y[t]) + \[Alpha]3 x[t] y[t] - c2 z[t] y[t];
eq3 = z'[t] == k1 - k2 z[t];

sol = NDSolveValue[{eq1, eq2, eq3, x[0] == x0, y[0] == y0, z[0] == z0}, {x, y, z}, {t, 0, 360}];

ParametricPlot3D[{sol[[1]][t], sol[[2]][t], sol[[3]][t]}, {t, 0, 360}, AxesLabel -> {x, y, z}, 
PlotRange -> All,  ImageSize -> {300, 300}],
{{c2, 0.4975, "c2"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny}, 
{{c1, 0.2215, "c1"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
{{k1, 0.8677, "k1"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny}, 
{{k2, 0.9611, "k2"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny}, 
{{\[Beta]2, 0.9817, "\[Beta]2"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny}, 
{{\[Alpha]3 , 0.2291, "\[Alpha]3"}, 0, 10, 0.1, Appearance -> "Labeled", ImageSize -> Tiny},
Delimiter,
{{x0, 1,"x[0]"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny}, 
{{y0, 1, "y[0]"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny}, 
{{z0, 0, "z[0]"}, 0, 10, .1, Appearance -> "Labeled", ImageSize -> Tiny}, 
                 ControlPlacement -> Left, ContinuousAction -> False]

enter image description here

Note: Experiment with different values for the parameters and the initial conditions.

$\endgroup$
  • $\begingroup$ THANKS. THAT IS GOOD FOR FIND A PARAMETRIC PLOT3D. WHAT ABOUT THE BIFRECTION AS SHOWN BY LINK@ZHK $\endgroup$ – sana alharbi Mar 7 at 4:58
  • $\begingroup$ @sanaalharbi You already got an answer for that by ChrisK $\endgroup$ – zhk Mar 7 at 5:04
2
$\begingroup$

I don't have time for a full analysis of your problem, but here are some thoughts and a rough approach to a bifurcation diagram.

First, the $V$ equation is decoupled from the others, so you should be able to solve for its equilibrium value $\hat V=k_1/k_2$ and inject it into the other two equations. Then regrouping of terms would reduce the number of independent parameters and also allow phase-plane analysis.

Once you've done that, I think you'll find the normal ($N_1$) and tumor cells ($T$) follow a combination of the Lotka-Volterra predator-prey model with logistic growth of both prey and predators. I think you could show that if a positive equilibrium exists it is unique and asymptotically stable.

This justifies a rough approach to generating a bifurcation diagram, using FindRoot to find an equilibrium for a given parameter value, then looping over the parameter value using the previous answer as an initial guess.

N1res = Tres = Vres = {}; (* results to be accumulated here *)
{N1i, Ti, Vi} = {21000, 0.6, k1/k2}; (* first initial guess *)
Do[
  eq = FindRoot[sys3 /. {_'[t] -> 0, var_[t] -> var},
    (* remove [t] and set time derivatives equal to zero *)
    {N1, N1i}, {T, Ti}, {V, Vi}];
  (* store results *)
  AppendTo[N1res, {k1, N1 /. eq}];
  AppendTo[Tres, {k1, T /. eq}];
  AppendTo[Vres, {k1, V /. eq}];
  (* update initial guess *)
  {N1i, Ti, Vi} = {N1, T, V} /. eq
, {k1, 0.6, 10.0, 0.1}]
GraphicsRow[{
  ListPlot[N1res, AxesLabel -> {"k1", "N1"}],
  ListPlot[Tres, AxesLabel -> {"k1", "T"}],
  ListPlot[Vres, AxesLabel -> {"k1", "V"}]
  }, ImageSize -> Large]

Mathematica graphics

This could be improved in many ways, but hopefully it gets you started.

$\endgroup$
  • $\begingroup$ I analyzed the stability of the model it was stable when the cells succeed to live together. but there is no growth for the third equation.@chris k $\endgroup$ – sana alharbi Mar 4 at 11:26
  • $\begingroup$ The equation for $V$ is independent of the other variables $N_1$ and $T$, so in all cases $V\to k_1/k_2$ as $t\to\infty$. $\endgroup$ – Chris K Mar 4 at 11:33
  • $\begingroup$ yes, that means, V is bounded, I proved that anlytically $\endgroup$ – sana alharbi Mar 4 at 15:32

Not the answer you're looking for? Browse other questions tagged or ask your own question.