How to Make a Bifurcation plot of a non linear system of ODEs

Question

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

Answer 1

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

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}]