1
$\begingroup$

I am currently working on a model of competition for resources between different species. I have 6 species competing for 3 resources. The first 3 species start at =0. the 4th species is added at t=1000 days, the 5th species is added at t=2000 days and the 6th species at t=5000. The first part is going well, but I can't seem to work out how I can get the last 3 species to be added later in the simulation. Can someone help me with this?

My code:

    (* Basic competition model *)
  (* Ns[i] \[Rule] population abundance of species i *)
  (* R[j] \[Rule] availability of resource j *)

  (* constant equation stuff: *)
  (* growth function of species i considering R of each resource at \
 time t *)
 \[Mu][i_, t_] := 
   Min[(r[[#]]*R[#, t])/(Ks[[#, i]] + R[#, t]) & /@ Range[k]];

 (* differential equations for Ns and R *)
 dNs[i_] := D[Ns[i, t], t] == Ns[i, t]*(\[Mu][i, t] - m[[i]]);
dR[j_] := 
  D[R[j, t], t] == 
   d*(S[[j]] - R[j, t]) - 
    Sum[c[[j, i]]*\[Mu][i, t]*Ns[i, t], {i, 1, n}];
(* Initialize parameters: 5 species, 5 resources, chaos *)
n = 6; (* number of species *)
k = 3; (* number of resources *)
tMax = 10000; (* time steps (days) *)
r = ConstantArray[1, n]; (* maximum growth rate of species *)
m = ConstantArray[0.25, n]; (* mortality rate of species *)
S = {6, 10, 14}; (* supply concentration of resources *)
d = 0.25; (* system turnover rate *)
Ks = ( {
   {1, 0.9, 0.3, 1.04, 0.34, 0.77},
   {0.3, 1, 0.9, 0.71, 1.02, 0.76},
   {0.9, 0.3, 1, 0.46, 0.34, 1.07}
   });  (* Subscript[K, ji ]= half saturation constant of resource j \
for species i *)
c = ({
   {0.04, 0.07, 0.04, 0.1, 0.03, 0.02},
   {0.08, 0.08, 0.1, 0.1, 0.05, 0.17},
   {0.14, 0.1, 0.1, 0.16, 0.06, 0.14}
    });  (* Subscript[c, ji] = content of resource j in species i *)
(* Ns and Rs at time t = 0 *)
initNs = (Ns[#, 0] == 0.1 + #/100) & /@ Range[n];

initR = (R[#, 0] == S[[#]]) & /@ Range[k];

(* make a list of n many Ns equations, and k many R equations *)
sysNs = Array[dNs[#] &, n];
sysR = Array[dR[#] &, k];
(* make lists of variables to be solved for *)
resultsNs = Array[Ns[#, t] &, n];
resultsRs = Array[R[#, t] &, n];
sol2 = NDSolveValue[{sysNs, sysR, initNs, initR}, {resultsNs, 
    resultsRs}, {t, 0, tMax}];
solNs2 = sol2[[1]];
solR2 = sol2[[2]];
Plot[solNs2, {t, 0, 10000}]        
$\endgroup$
  • 1
    $\begingroup$ Would WhenEvent help? sol2 = NDSolveValue[{sysNs, sysR, initNs, initR, WhenEvent[t > 1000, Ns[4, t] -> 100]}, {resultsNs, resultsRs}, {t, 0, tMax}] $\endgroup$ – Ruud3.1415 Oct 2 at 10:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.