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

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

• Could you explain the goal of the modification? As I see, if you model two independent logistic equation there will always look as the panel on the right. The interaction between the populations in the first set of equations is given by dS/dt which depends of both populations, while each populations only depends on S and it self. Jun 16, 2021 at 19:24
• The point is: the implicit form is present in more examples than the explicit, thus it would be the preferred way to analyze the interactions. By converting the first set into the more common way of representation, I will have a consistent data set instead of having some described in one way and other in another... Jun 17, 2021 at 6:13
• I would show an example of a numerical solution, but you don't give values for D and S0 ... Jun 17, 2021 at 14:12
• Thanks, I have updated the question to incorporate these values. The example is taken from SCIENCE 1980;207:1491 (E. coli C-8 = N₁; P. aeruginosa PA0283 = N₂). Jun 18, 2021 at 13:07
• in DOI: 10.1126/science.192.4238.463 I have found the relation μ = μmax[S/(S + K)]. This allows to convert d/dt N₁ = (μ₁SN₁)/(K₁+S) to d/dt N₁ = μN₁. Only the logistic term is missing, which I can guess estimating the carrying capacity. Jun 18, 2021 at 15:46

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()

• Probably, one could also consider the limit where the neutrient dynamics is very slow? Jun 18, 2021 at 13:42
• Thanks, I'll need some time to digest this answer... Jun 19, 2021 at 15:16