2
$\begingroup$

I am modelling bacterial growth. I have a quite difficult non linear systems of ODE's. I have 2 bacteria that i'm modelling with competition, as wel as transfer to a different place. so there are 2 ODEs for every bacteria because when it moves it will be a different one (P en S equation). Now i am trying to do a parameter analysis and i've been trying for days to make a bifurcation plot with on the x-axis one of the parameters (a0,b0,d1,k,awmin,tmin). But nothing i've been trying so far worked. Can you guys help me out a little and show how it works? here is the code of the bacteria, the bifurcation should take a spot around t=150, I've had a lot of trouble with this.

sp1 = D[G[t], t] == 
   Growth*G[t]*((Sp - G[t] - a0*B[t] )/Sp) + P - S  ;  
sp2 = D[M[t], 
    t] == (Growth2 - Kill)*M[t]*((Ss - M[t] - a0*A[t] )/Ss) + S - 
    P   ; 
P := M[t]* d1 ;  
S := G[t]*s*((G[t] + a0*B[t])/Sp) + G[t]*k     ;  
sp3 = D[B[t], t] == Growth*B[t]*((Sp - B[t] - b0*G[t] )/Sp) + P2 - S2;
sp4 = D[A[t], t] == 
   Growth2* A[t]*((Ss - A[t] - b0*M[t] )/Ss) + S2 - P2;
P2 := A[t]* d1;       
S2 := B[t]*0.01*((B[t] + b0*G[t] )/Sp) + B[t]*k;    
(*Parameters*)

Growth = 14.8*((T - tmin)*(Aw - Amin))^2  ;  
Growth2 = 7.4*((T - tmin)*(Aw - Amin))^2  ;
T := 4 + 18.09*Sin[0.01016*t + 0.3418];
Aw := 0.97 + 0.0001*t   ;
Kill = 1.5;
tmin = 12;
Amin = 0.95;
Sp = 1000;
Ss = 500;
a0 = 0.4;
b0 = 0.21;
d1 = 0.05;
k = 0.01;
s = 0.001;
all = { sp1, sp2, sp3, sp4};  
init1 = {G[0] == 60, M[0] == 3, B[0] == 50, A[0] == 30}; 
Solution = 
  NDSolveValue[{all, init1}, {G[t], M[t], B[t], A[t]}, {t, 0, 200}];  
Plot[Solution, {t, 0, 200}, PlotStyle -> {Red, Pink, Blue, Cyan}] 
ParametricPlot[{Solution[[1]], Solution[[3]]}, {t, 0, 200}, 
 PlotRange -> {{0, 100}, {0, 100}}, PlotStyle -> Red, 
 AspectRatio -> 1, PlotLabel -> "Phase plot", AxesLabel -> {"G", "B"}]

Plot and parametric plot of the ODE's

$\endgroup$
1
  • 3
    $\begingroup$ This doesn't fix your problem, but you'd be better off if you rewrote your code in a more ordered fashion like this pastebin.com/LWVXcmWC where definitions come before usage and things are segregated into functions of t vs. parameters. $\endgroup$
    – flinty
    Aug 21, 2020 at 17:42

1 Answer 1

2
$\begingroup$

We can use Module for parametric research as follows

f[a0p_, b0p_, d1p_, kp_, Ap_, tp_] := 
 Module[{a0 = a0p, b0 = b0p, d1 = d1p, k = kp, Amin = Ap, tmin = tp}, 
  Kill = 1.5;
  (*tmin=12;
  Amin=0.95;*)
  Sp = 1000;
  Ss = 500;
  (*a0=0.4;
  b0=0.21;
  d1=0.05;
  k=0.01;*)
  s = 0.001; 
  sp1 = D[G[t], t] == Growth*G[t]*((Sp - G[t] - a0*B[t])/Sp) + P - S;
  sp2 = D[M[t], 
     t] == (Growth2 - Kill)*M[t]*((Ss - M[t] - a0*A[t])/Ss) + S - P;
  P := M[t]*d1;
  S := G[t]*s*((G[t] + a0*B[t])/Sp) + G[t]*k;
  sp3 = D[B[t], t] == Growth*B[t]*((Sp - B[t] - b0*G[t])/Sp) + P2 - S2;
  sp4 = D[A[t], t] == 
    Growth2*A[t]*((Ss - A[t] - b0*M[t])/Ss) + S2 - P2;
  P2 := A[t]*d1;
  S2 := B[t]*0.01*((B[t] + b0*G[t])/Sp) + B[t]*k; 
  Growth = 14.8*((T - tmin)*(Aw - Amin))^2;
  Growth2 = 7.4*((T - tmin)*(Aw - Amin))^2;
  T := 4 + 18.09*Sin[0.01016*t + 0.3418];
  Aw := 0.97 + 0.0001*t; all = {sp1, sp2, sp3, sp4};
  init1 = {G[0] == 60, M[0] == 3, B[0] == 50, A[0] == 30};
  Solution = 
   NDSolveValue[{all, init1}, {G[t], M[t], B[t], A[t]}, {t, 0, 200}]; 
  Solution]

To plot Solution we use f as

sol = f[.4, .21, .05, .01, .95, 12];
{Plot[sol, {t, 0, 200}, PlotStyle -> {Red, Pink, Blue, Cyan}],
 ParametricPlot[{sol[[1]], sol[[3]]}, {t, 0, 200}, 
  PlotRange -> {{0, 100}, {0, 100}}, PlotStyle -> Red, 
  AspectRatio -> 1, PlotLabel -> "Phase plot", 
  AxesLabel -> {"G", "B"}]}

Figure1

To plot Solution as a function of parameters for a given time t=150 we define new function, for example,

solp[ap_] := f[ap, .21, .05, .01, .95, 12] /. t -> 150;
solp1[bp_] := f[.4, bp, .05, .01, .95, 12] /. t -> 150;
var = {G, M, B, A}; Table[
 Plot[solp[ap][[i]], {ap, .1, .5}, Frame -> True, 
  FrameLabel -> {"a0", var[[i]]}], {i, 4}]


Table[Plot[solp1[bp][[i]], {bp, .1, .5}, Frame -> True, 
  FrameLabel -> {"b0", var[[i]]}], {i, 4}]

Figure 1

$\endgroup$
4
  • 1
    $\begingroup$ I think ParametricNDSolve(Value)[] might be more conducive to exploring parameter spaces. $\endgroup$ Aug 23, 2020 at 14:50
  • $\begingroup$ @J.M. Yes, we can use ParametricNDSolve. It gives actually the same result. $\endgroup$ Aug 23, 2020 at 15:06
  • $\begingroup$ Thank you so much for your response, it is very useful! You just made my day :D $\endgroup$
    – Zara
    Aug 24, 2020 at 12:31
  • $\begingroup$ @Zara You are welcome! $\endgroup$ Aug 24, 2020 at 14:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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