How to convert system of equations with explicit substrate consumption into an implicit one?

Question

I have a description of competition between two bacteria that is described by the system of equations:

S = substrate, μ = bacterial division rate, K = half saturation constant, y = yield (cells per unit of nutrient), N = concentration of bacteria.

Is it possible to convert it to a system of implicit substrate consumption such as the logistic system? For instance into this:

where μ = maximum growth rate and K = carrying capacity.

Essentially, the term Ṡ should disappear, μⱼS become μ, and Kₛ+S become K. I tried such conversion but I got the figure on the right panel instead of the original on the left:

On the left,

μ = μ₁ = 0.81
ν = μ₂ = 0.91
y₁ = 2.5*10^10 
y₂ = 3.8*10^10
K₁ = 3.0*10^-6 
K₂ = 3.1*10^-4
S₀ = 1*10^-4 (g/L)
N¹₀ = 1.0*10^2    
N²₀ = N¹₀*200     
D = 6*10^-2 (/h)
Vol = 200 mL

Obviously, the conversion I did has not included the terms y₁, y₂, K₁, and K₂. Thus, the red species never took over. How am I to bind all these variables in fewer terms into a logistic system?

Thank you

Answer 1

In general, no. You're trying to reduce the a three-dimensional dynamical system ({S, N1, N2}) to a two-dimensional system ({N1, N2}), which can't be done in general.

One way to approximate a two-dimensional system is to assume that the nutrient dynamics are on a fast time scale, so that for a given N1, N2 you can write down the quasi-equilibrium S*(N1,N2) and plug it back into the N1 and N2 equations. If I've done the math right, the value of S* as a function of N1, N2, and the other model parameters is the solution to the cubic equation dS/dt = 0, or:

for S. Unfortunately, this probably means the answer (although it will have a closed-form solution) will be too complex to do anything with analytically, although you could solve it computationally.

I tried this (code below): the dynamics actually track the full system very well for a while, but then the equations seem to blow up (the lines are the full 3-d dynamics, the points are the quasi-equilibrium approximation). It would take a lot more thought and analysis to figure out what's going on here — most likely, the equations just aren't biologically well-posed. (In the full system, S crashes down from about 1e-4 to 2e-7 very suddenly around t=12, then declines extremely slowly after that ...)

@RogerVadim suggested that one could also try the limit of slow resource dynamics, but in this case (assuming we mean dS/dt << 1) the dynamics are simple and unrealistic; slow resource dynamics means we would hold S constant at its starting value and see what happened. In this case, because the per-capita growth rates of N1 and N2 are constant at constant S, we would just get exponential growth or decay for each species (depending on whether the birth rate was greater or less than the death rate). There are plenty of other things you could try to simplify the system - what if D=0? What if D is non-zero but << 1 ? Is there a way fix the 2-D (fast resource dynamics) system so that it works, either by changing something computational or reformulating the system? - but nothing easy is immediately obvious to me.

p0 <- c(mu1 = 0.81, mu2 = 0.91, y1 = 2.5e10,
        y2 = 3.8e10, K1 = 3.0e-6, K2 = 3.1e-4,
        S0 = 1e-4, D=6e-2)

y0 <- c(S = 1e-4, N1 = 1e2, N2 = 2e4)

library(deSolve)
g1 <- function(t,y,parms) {
  g <- with(as.list(c(y,parms)),
  {
    m1 <- mu1/y1*(S*N1)/(K1 + S)
    m2 <- mu2/y2*(S*N2)/(K2 + S)g
    c(S = (S0-S)*D-m1-m2,
      N1 = m1*y1-D*N1,
      N2 = m2*y2-D*N2)
  }
  )
  list(g)
}

Seqfun <- function(S, y, parms) {
  Seq <- with(as.list(c(y, parms)),
  (S0-S)*D - mu1/y1*N1*S/(K1+S) - mu2/y2*N2*S/(K2+S))
  return(Seq)
}


g2 <- function(t,y,parms) {
  Seq <- uniroot(Seqfun, interval = c(0, 1e3), y = y, parms = parms)$root
  cat(y[["N1"]], y[["N2"]], Seq, "\n")
  g <- with(as.list(c(y, parms)),
            c(N1 = mu1*(Seq*N1)/(K1 + Seq) - D*N1,
              N2 = mu2*(Seq*N2)/(K2 + Seq) - D*N2)
            )
  list(g, S = Seq)
}

tvec <- seq(0, 60)
r1 <- ode(y = y0, times = tvec, func = g1, parms = p0)
r2 <- ode(y = y0[-1], times = tvec, func = g2, parms = p0)

png("bact.png")
par(las=1, bty="l")
matplot(r1[,-(1:2)], type="l", log="y", lty =1, xlab="time", ylab="density")
matpoints(r2[,2:3], pch=1)
dev.off()