from sympy import *
init_printing()

"""
Explanation: R1C1 Ecosystems Model Equilibria
In this notebook we explore the equilibria of the one resource / one consumer (R1-C1) ecosystem model.
End of explanation
"""


x = Symbol('x') # biomass of resouce species [biomassX]
y = Symbol('y') # biomass of consumer species [biomassY]

"""
Explanation: We start by defining the variables:
End of explanation
"""


r = Symbol('r') # growth rate [biomassX/(biomassX*second)]
k = Symbol('k') # carrying capacity [biomassX]
a = Symbol('a') # attack rate [biomassY/biomassX]
d = Symbol('d') # death rate due to predation [biomassX/(biomassY*second)]
e = Symbol('e') # conversion efficiency [biomassY/(biomassY*second)]
h = Symbol('h') # half saturaion [biomassY]
m = Symbol('m') # mortality [biomassY/(biomassY*second)]

"""
Explanation: We now define the parameters:
End of explanation
"""


deltaX = r * x * (1 - (x / k)) - d * (a * y * x)/(h+a*x) # [biomassX/second]
deltaY = (e * a * y * x)/(h+a*x) - m * y                 # [biomassY/second]
deltaX, deltaY

"""
Explanation: Finally, we define the r1-c1 equations:
End of explanation
"""


equilibria = solve([deltaX, deltaY], [x, y])
equilibria

"""
Explanation: Solving these equations for the equilibria we get:
End of explanation
"""


cX0 = equilibria[2][0]
cY0 = equilibria[2][1]
(cX0, cY0)

factor(deltaX.subs({ x: cX0, y: cY0}))

factor(deltaY.subs({ x: cX0, y: cY0}))

"""
Explanation: The ecosystem with no resources and hence no consumers, as well as the ecosystem with only resources at their carrying capacity, are the first two, fairly obvious equilibria.
The more interesting equilibria is the third one which is highly likely to be the result of a hopf-bifurcation and the focus of the cyclycial dynamics of the overall ecosystem.
We define, $cX$ and $cY$ to be:
End of explanation
"""


hHat = Symbol('\hat{h}')
hHatDef = h/(e-m)
# [biomassY*second]
#   == [biomassY]/([biomassY/(biomassY*second)] - [biomassY/(biomassY*second)])
#   == [biomassY]/([biomassY/(biomassY*second)])
#   == [biomassY] * [(biomassY*second)/biomassY]  
hHat, hHatDef

"""
Explanation: In these equations, the collection of the parameters, $h, a$, and $(e-m)$ seems to be an important value. To explore this, we define the constant, $\hat{h}$, as a scaled half-saturation, as follows:
End of explanation
"""


cX1 = cX0.subs(hHatDef, hHat) # [biomassX] == [seconds/biomassX]*[biomassY/(biomassY*second)]
cX1

"""
Explanation: Applying this simplification to $cX$ we get
End of explanation
"""


cY0.subs(hHatDef, hHat)

"""
Explanation: Unfortunately, with out some additional work, we can only partially apply this simplificaiton to $cY$. An initial partial application gives us:
End of explanation
"""


cYfactor0 = (-a*e*k + a*k*m + h*m)/(a*k*(e-m))
cYfactor0

"""
Explanation: To simplify this further we need to focus upon the factor:
End of explanation
"""


cYfactor1 = h*m/(a*k*(e-m)) - 1
cYfactor1

"""
Explanation: We assert that $cYfactor0$ can be rewritten as:
End of explanation
"""


expand(cYfactor0)

expand(factor(cYfactor1))

"""
Explanation: To prove this we note that:
End of explanation
"""


expand(cYfactor0) == expand(factor(cYfactor1))

"""
Explanation: Hence:
End of explanation
"""


cY0.subs(hHatDef, hHat).subs(cYfactor0, cYfactor1).subs(hHatDef, hHat)

"""
Explanation: Applying this substitution to $cY$ we get:
End of explanation
"""


cY1 = ((hHat*e*r)/(a*d))*(1-hHat*m/(a*k))
cY1

"""
Explanation: or after multiplying the $-1$ into $(\frac{\hat{h}m}{ak} - 1)$, we obtain:
End of explanation
"""


cX1, cY1

"""
Explanation: Giving us the resulting cyclic equilibria:
End of explanation
"""


deltaX, deltaY

"""
Explanation: We note that this is only an ecologically relevant equilibria if both $0 < \hat{h}$ and $\frac{\hat{h}m}{ak} < 1$. Given the definition of $\hat{h} = \frac{h}{e-m}$, we know that $0 < \hat{h}$, only happens when $m < e$. Similarly, $\frac{\hat{h}m}{ak} < 1$ only when $\hat{h} < \frac{ak}{m}$.
Stability about each equilibria
We now study the stability of each equlibria by looking at the jacobian of the r1-c1 equations. Recall that the equations themselves are:
End of explanation
"""


diff(deltaX, x), diff(deltaX, y)

diff(deltaY,x), diff(deltaY,y)

from sympy.matrices import *
delta = Matrix([deltaX, deltaY])
jacob = delta.jacobian([x, y])
jacob

"""
Explanation: The jacobian is then:
End of explanation
"""


jacob.subs(x, 0).subs(y, 0)

jacob.subs(x, 0).subs(y, 0).det()

jacob.subs(x, 0).subs(y, 0).eigenvects()

"""
Explanation: Stability around origin (no ecology)
End of explanation
"""


jacob.subs(x, k).subs(y, 0)

"""
Explanation: Stability around resources at carrying capacity
End of explanation
"""


aHat = Symbol('\hat{a}')
aHatDef = a*k/(a*k + h)
#  [??] == ??
aHat, aHatDef

jacob.subs(x, k).subs(y, 0).subs(aHatDef, aHat).det()

jacob.subs(x, k).subs(y, 0).subs(aHatDef, aHat).eigenvects()

"""
Explanation: We consider a scaled attack rate:
End of explanation
"""


jacob.subs(x, cX1).subs(y, cY1)

"""
Explanation: Stability around cyclic equilibrium
End of explanation
"""


bHat = Symbol('\hat{b}')
bHatDef = 1 - hHat*m/(a*k)
#  [??] == ??
bHat, bHatDef

"""
Explanation: Consider a scaled ??:
End of explanation
"""


cHat = Symbol('\hat{c}')
cHatDef = hHat/(hHat*m+h)
#  [??] == ??
cHat, cHatDef

jacob.subs(x, cX1).subs(y, cY1).subs(bHatDef, bHat).subs(cHatDef, cHat)

jacob.subs(x, cX1).subs(y, cY1).subs(bHatDef, bHat).det()

jacob.subs(x, cX1).subs(y, cY1).eigenvects()

"""
Explanation: and a scaled ??:
End of explanation
"""


