R1C1 Ecosystems Model Equilibria

In this notebook we explore the equilibria of the one resource / one consumer (R1-C1) ecosystem model.



In [1]:

    
from sympy import *
init_printing()

We start by defining the variables:



In [2]:

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

We now define the parameters:



In [3]:

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

Finally, we define the r1-c1 equations:



In [4]:

    
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









    Out[4]:





$$\left ( - \frac{a d x y}{a x + h} + r x \left(1 - \frac{x}{k}\right), \quad \frac{a e x y}{a x + h} - m y\right )$$

Solving these equations for the equilibria we get:



In [5]:

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









    Out[5]:





$$\left [ \left ( 0, \quad 0\right ), \quad \left ( k, \quad 0\right ), \quad \left ( \frac{h m}{a \left(e - m\right)}, \quad - \frac{e h r}{a^{2} d k \left(e - m\right)^{2}} \left(- a e k + a k m + h m\right)\right )\right ]$$

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:



In [6]:

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









    Out[6]:





$$\left ( \frac{h m}{a \left(e - m\right)}, \quad - \frac{e h r}{a^{2} d k \left(e - m\right)^{2}} \left(- a e k + a k m + h m\right)\right )$$



In [7]:

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









    Out[7]:





$$0$$



In [8]:

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









    Out[8]:





$$0$$

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:



In [9]:

    
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









    Out[9]:





$$\left ( \hat{h}, \quad \frac{h}{e - m}\right )$$

Applying this simplification to $cX$ we get



In [10]:

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









    Out[10]:





$$\frac{\hat{h} m}{a}$$

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



In [11]:

    
cY0.subs(hHatDef, hHat)









    Out[11]:





$$- \frac{\hat{h} e r}{a^{2} d k \left(e - m\right)} \left(- a e k + a k m + h m\right)$$

To simplify this further we need to focus upon the factor:



In [12]:

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









    Out[12]:





$$\frac{1}{a k \left(e - m\right)} \left(- a e k + a k m + h m\right)$$

We assert that $cYfactor0$ can be rewritten as:



In [13]:

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









    Out[13]:





$$-1 + \frac{h m}{a k \left(e - m\right)}$$

To prove this we note that:



In [14]:

    
expand(cYfactor0)









    Out[14]:





$$- \frac{a e k}{a e k - a k m} + \frac{a k m}{a e k - a k m} + \frac{h m}{a e k - a k m}$$



In [15]:

    
expand(factor(cYfactor1))









    Out[15]:





$$- \frac{a e k}{a e k - a k m} + \frac{a k m}{a e k - a k m} + \frac{h m}{a e k - a k m}$$

Hence:



In [16]:

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









    Out[16]:





True

Applying this substitution to $cY$ we get:



In [17]:

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









    Out[17]:





$$- \frac{\hat{h} e r}{a d} \left(\frac{\hat{h} m}{a k} - 1\right)$$

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



In [18]:

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









    Out[18]:





$$\frac{\hat{h} e r}{a d} \left(- \frac{\hat{h} m}{a k} + 1\right)$$

Giving us the resulting cyclic equilibria:



In [19]:

    
cX1, cY1









    Out[19]:





$$\left ( \frac{\hat{h} m}{a}, \quad \frac{\hat{h} e r}{a d} \left(- \frac{\hat{h} m}{a k} + 1\right)\right )$$

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:



In [20]:

    
deltaX, deltaY









    Out[20]:





$$\left ( - \frac{a d x y}{a x + h} + r x \left(1 - \frac{x}{k}\right), \quad \frac{a e x y}{a x + h} - m y\right )$$

The jacobian is then:



In [21]:

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









    Out[21]:





$$\left ( \frac{a^{2} d x y}{\left(a x + h\right)^{2}} - \frac{a d y}{a x + h} + r \left(1 - \frac{x}{k}\right) - \frac{r x}{k}, \quad - \frac{a d x}{a x + h}\right )$$



In [22]:

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









    Out[22]:





$$\left ( - \frac{a^{2} e x y}{\left(a x + h\right)^{2}} + \frac{a e y}{a x + h}, \quad \frac{a e x}{a x + h} - m\right )$$



In [23]:

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









    Out[23]:





$$\left[\begin{matrix}\frac{a^{2} d x y}{\left(a x + h\right)^{2}} - \frac{a d y}{a x + h} + r \left(1 - \frac{x}{k}\right) - \frac{r x}{k} & - \frac{a d x}{a x + h}\\- \frac{a^{2} e x y}{\left(a x + h\right)^{2}} + \frac{a e y}{a x + h} & \frac{a e x}{a x + h} - m\end{matrix}\right]$$

Stability around origin (no ecology)



In [24]:

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









    Out[24]:





$$\left[\begin{matrix}r & 0\\0 & - m\end{matrix}\right]$$



In [25]:

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









    Out[25]:





$$- m r$$



In [26]:

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









    Out[26]:





$$\left [ \left ( - m, \quad 1, \quad \left [ \left[\begin{matrix}0\\1\end{matrix}\right]\right ]\right ), \quad \left ( r, \quad 1, \quad \left [ \left[\begin{matrix}1\\0\end{matrix}\right]\right ]\right )\right ]$$

Stability around resources at carrying capacity



In [27]:

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









    Out[27]:





$$\left[\begin{matrix}- r & - \frac{a d k}{a k + h}\\0 & \frac{a e k}{a k + h} - m\end{matrix}\right]$$

We consider a scaled attack rate:



In [28]:

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









    Out[28]:





$$\left ( \hat{a}, \quad \frac{a k}{a k + h}\right )$$



In [29]:

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









    Out[29]:





$$- \hat{a} e r + m r$$



In [30]:

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









    Out[30]:





$$\left [ \left ( - r, \quad 1, \quad \left [ \left[\begin{matrix}1\\0\end{matrix}\right]\right ]\right ), \quad \left ( \hat{a} e - m, \quad 1, \quad \left [ \left[\begin{matrix}\frac{\hat{a} d}{- \hat{a} e + m - r}\\1\end{matrix}\right]\right ]\right )\right ]$$

Stability around cyclic equilibrium



In [31]:

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









    Out[31]:





$$\left[\begin{matrix}\frac{\hat{h}^{2} e m r}{\left(\hat{h} m + h\right)^{2}} \left(- \frac{\hat{h} m}{a k} + 1\right) - \frac{\hat{h} e r}{\hat{h} m + h} \left(- \frac{\hat{h} m}{a k} + 1\right) - \frac{\hat{h} m r}{a k} + r \left(- \frac{\hat{h} m}{a k} + 1\right) & - \frac{\hat{h} d m}{\hat{h} m + h}\\- \frac{\hat{h}^{2} e^{2} m r \left(- \frac{\hat{h} m}{a k} + 1\right)}{d \left(\hat{h} m + h\right)^{2}} + \frac{\hat{h} e^{2} r \left(- \frac{\hat{h} m}{a k} + 1\right)}{d \left(\hat{h} m + h\right)} & \frac{\hat{h} e m}{\hat{h} m + h} - m\end{matrix}\right]$$

Consider a scaled ??:



In [32]:

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









    Out[32]:





$$\left ( \hat{b}, \quad - \frac{\hat{h} m}{a k} + 1\right )$$

and a scaled ??:



In [33]:

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









    Out[33]:





$$\left ( \hat{c}, \quad \frac{\hat{h}}{\hat{h} m + h}\right )$$



In [34]:

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









    Out[34]:





$$\left[\begin{matrix}\hat{b} \hat{c}^{2} e m r - \hat{b} \hat{c} e r + \hat{b} r - \frac{\hat{h} m r}{a k} & - \hat{c} d m\\- \frac{\hat{b} m}{d} \hat{c}^{2} e^{2} r + \frac{\hat{b} \hat{c}}{d} e^{2} r & \hat{c} e m - m\end{matrix}\right]$$



In [35]:

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









    Out[35]:





$$- \frac{\hat{b} \hat{h}^{3} d e^{2} m^{2} r}{\hat{h}^{3} d m^{3} + 3 \hat{h}^{2} d h m^{2} + 3 \hat{h} d h^{2} m + d h^{3}} + \frac{\hat{b} \hat{h}^{3} e^{2} m^{2} r}{\hat{h}^{3} m^{3} + 3 \hat{h}^{2} h m^{2} + 3 \hat{h} h^{2} m + h^{3}} + \frac{\hat{b} \hat{h}^{2} d e^{2} m r}{\hat{h}^{2} d m^{2} + 2 \hat{h} d h m + d h^{2}} - \frac{\hat{b} \hat{h}^{2} e^{2} m r}{\hat{h}^{2} m^{2} + 2 \hat{h} h m + h^{2}} - \frac{\hat{b} \hat{h}^{2} e m^{2} r}{\hat{h}^{2} m^{2} + 2 \hat{h} h m + h^{2}} + \frac{2 \hat{b} \hat{h} e m r}{\hat{h} m + h} - \hat{b} m r - \frac{\hat{h}^{2} e m^{2} r}{\hat{h} a k m + a h k} + \frac{\hat{h} m^{2} r}{a k}$$



In [36]:

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









    Out[36]:





$$\left [ \left ( - \frac{1}{2 a k \left(\hat{h}^{2} m^{2} + 2 \hat{h} h m + h^{2}\right)} \sqrt{4 \hat{h}^{6} m^{6} r^{2} + 4 \hat{h}^{5} a e k m^{5} r - 4 \hat{h}^{5} a k m^{6} r - 4 \hat{h}^{5} a k m^{5} r^{2} - 4 \hat{h}^{5} e h m^{4} r^{2} + 16 \hat{h}^{5} h m^{5} r^{2} + \hat{h}^{4} a^{2} e^{2} k^{2} m^{4} - 2 \hat{h}^{4} a^{2} e k^{2} m^{5} - 2 \hat{h}^{4} a^{2} e k^{2} m^{4} r + \hat{h}^{4} a^{2} k^{2} m^{6} + 2 \hat{h}^{4} a^{2} k^{2} m^{5} r + \hat{h}^{4} a^{2} k^{2} m^{4} r^{2} + 2 \hat{h}^{4} a e^{2} h k m^{3} r + 14 \hat{h}^{4} a e h k m^{4} r + 6 \hat{h}^{4} a e h k m^{3} r^{2} - 16 \hat{h}^{4} a h k m^{5} r - 16 \hat{h}^{4} a h k m^{4} r^{2} + \hat{h}^{4} e^{2} h^{2} m^{2} r^{2} - 8 \hat{h}^{4} e h^{2} m^{3} r^{2} + 24 \hat{h}^{4} h^{2} m^{4} r^{2} + 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{3} - 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{2} r - 6 \hat{h}^{3} a^{2} e h k^{2} m^{4} - 8 \hat{h}^{3} a^{2} e h k^{2} m^{3} r - 2 \hat{h}^{3} a^{2} e h k^{2} m^{2} r^{2} + 4 \hat{h}^{3} a^{2} h k^{2} m^{5} + 8 \hat{h}^{3} a^{2} h k^{2} m^{4} r + 4 \hat{h}^{3} a^{2} h k^{2} m^{3} r^{2} + 2 \hat{h}^{3} a e^{2} h^{2} k m^{2} r - 2 \hat{h}^{3} a e^{2} h^{2} k m r^{2} + 16 \hat{h}^{3} a e h^{2} k m^{3} r + 12 \hat{h}^{3} a e h^{2} k m^{2} r^{2} - 24 \hat{h}^{3} a h^{2} k m^{4} r - 24 \hat{h}^{3} a h^{2} k m^{3} r^{2} - 4 \hat{h}^{3} e h^{3} m^{2} r^{2} + 16 \hat{h}^{3} h^{3} m^{3} r^{2} + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m^{2} - 2 \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m r + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} r^{2} - 6 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{3} - 10 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{2} r - 4 \hat{h}^{2} a^{2} e h^{2} k^{2} m r^{2} + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{4} + 12 \hat{h}^{2} a^{2} h^{2} k^{2} m^{3} r + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{2} r^{2} + 6 \hat{h}^{2} a e h^{3} k m^{2} r + 6 \hat{h}^{2} a e h^{3} k m r^{2} - 16 \hat{h}^{2} a h^{3} k m^{3} r - 16 \hat{h}^{2} a h^{3} k m^{2} r^{2} + 4 \hat{h}^{2} h^{4} m^{2} r^{2} - 2 \hat{h} a^{2} e h^{3} k^{2} m^{2} - 4 \hat{h} a^{2} e h^{3} k^{2} m r - 2 \hat{h} a^{2} e h^{3} k^{2} r^{2} + 4 \hat{h} a^{2} h^{3} k^{2} m^{3} + 8 \hat{h} a^{2} h^{3} k^{2} m^{2} r + 4 \hat{h} a^{2} h^{3} k^{2} m r^{2} - 4 \hat{h} a h^{4} k m^{2} r - 4 \hat{h} a h^{4} k m r^{2} + a^{2} h^{4} k^{2} m^{2} + 2 a^{2} h^{4} k^{2} m r + a^{2} h^{4} k^{2} r^{2}} - \frac{1}{2 a k \left(\hat{h} m + h\right)^{2}} \left(2 \hat{h}^{3} m^{3} r - \hat{h}^{2} a e k m^{2} + \hat{h}^{2} a k m^{3} - \hat{h}^{2} a k m^{2} r - \hat{h}^{2} e h m r + 4 \hat{h}^{2} h m^{2} r - \hat{h} a e h k m + \hat{h} a e h k r + 2 \hat{h} a h k m^{2} - 2 \hat{h} a h k m r + 2 \hat{h} h^{2} m r + a h^{2} k m - a h^{2} k r\right), \quad 1, \quad \left [ \left[\begin{matrix}\frac{2 \hat{h} a d k m \left(\hat{h}^{2} m^{2} + 2 \hat{h} h m + h^{2}\right)}{\left(\hat{h} m + h\right) \left(- 2 \hat{h}^{3} m^{3} r - \hat{h}^{2} a e k m^{2} + \hat{h}^{2} a k m^{3} + \hat{h}^{2} a k m^{2} r + \hat{h}^{2} e h m r - 4 \hat{h}^{2} h m^{2} r - \hat{h} a e h k m - \hat{h} a e h k r + 2 \hat{h} a h k m^{2} + 2 \hat{h} a h k m r - 2 \hat{h} h^{2} m r + a h^{2} k m + a h^{2} k r + \sqrt{4 \hat{h}^{6} m^{6} r^{2} + 4 \hat{h}^{5} a e k m^{5} r - 4 \hat{h}^{5} a k m^{6} r - 4 \hat{h}^{5} a k m^{5} r^{2} - 4 \hat{h}^{5} e h m^{4} r^{2} + 16 \hat{h}^{5} h m^{5} r^{2} + \hat{h}^{4} a^{2} e^{2} k^{2} m^{4} - 2 \hat{h}^{4} a^{2} e k^{2} m^{5} - 2 \hat{h}^{4} a^{2} e k^{2} m^{4} r + \hat{h}^{4} a^{2} k^{2} m^{6} + 2 \hat{h}^{4} a^{2} k^{2} m^{5} r + \hat{h}^{4} a^{2} k^{2} m^{4} r^{2} + 2 \hat{h}^{4} a e^{2} h k m^{3} r + 14 \hat{h}^{4} a e h k m^{4} r + 6 \hat{h}^{4} a e h k m^{3} r^{2} - 16 \hat{h}^{4} a h k m^{5} r - 16 \hat{h}^{4} a h k m^{4} r^{2} + \hat{h}^{4} e^{2} h^{2} m^{2} r^{2} - 8 \hat{h}^{4} e h^{2} m^{3} r^{2} + 24 \hat{h}^{4} h^{2} m^{4} r^{2} + 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{3} - 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{2} r - 6 \hat{h}^{3} a^{2} e h k^{2} m^{4} - 8 \hat{h}^{3} a^{2} e h k^{2} m^{3} r - 2 \hat{h}^{3} a^{2} e h k^{2} m^{2} r^{2} + 4 \hat{h}^{3} a^{2} h k^{2} m^{5} + 8 \hat{h}^{3} a^{2} h k^{2} m^{4} r + 4 \hat{h}^{3} a^{2} h k^{2} m^{3} r^{2} + 2 \hat{h}^{3} a e^{2} h^{2} k m^{2} r - 2 \hat{h}^{3} a e^{2} h^{2} k m r^{2} + 16 \hat{h}^{3} a e h^{2} k m^{3} r + 12 \hat{h}^{3} a e h^{2} k m^{2} r^{2} - 24 \hat{h}^{3} a h^{2} k m^{4} r - 24 \hat{h}^{3} a h^{2} k m^{3} r^{2} - 4 \hat{h}^{3} e h^{3} m^{2} r^{2} + 16 \hat{h}^{3} h^{3} m^{3} r^{2} + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m^{2} - 2 \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m r + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} r^{2} - 6 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{3} - 10 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{2} r - 4 \hat{h}^{2} a^{2} e h^{2} k^{2} m r^{2} + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{4} + 12 \hat{h}^{2} a^{2} h^{2} k^{2} m^{3} r + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{2} r^{2} + 6 \hat{h}^{2} a e h^{3} k m^{2} r + 6 \hat{h}^{2} a e h^{3} k m r^{2} - 16 \hat{h}^{2} a h^{3} k m^{3} r - 16 \hat{h}^{2} a h^{3} k m^{2} r^{2} + 4 \hat{h}^{2} h^{4} m^{2} r^{2} - 2 \hat{h} a^{2} e h^{3} k^{2} m^{2} - 4 \hat{h} a^{2} e h^{3} k^{2} m r - 2 \hat{h} a^{2} e h^{3} k^{2} r^{2} + 4 \hat{h} a^{2} h^{3} k^{2} m^{3} + 8 \hat{h} a^{2} h^{3} k^{2} m^{2} r + 4 \hat{h} a^{2} h^{3} k^{2} m r^{2} - 4 \hat{h} a h^{4} k m^{2} r - 4 \hat{h} a h^{4} k m r^{2} + a^{2} h^{4} k^{2} m^{2} + 2 a^{2} h^{4} k^{2} m r + a^{2} h^{4} k^{2} r^{2}}\right)}\\1\end{matrix}\right]\right ]\right ), \quad \left ( \frac{1}{2 a k \left(\hat{h}^{2} m^{2} + 2 \hat{h} h m + h^{2}\right)} \sqrt{4 \hat{h}^{6} m^{6} r^{2} + 4 \hat{h}^{5} a e k m^{5} r - 4 \hat{h}^{5} a k m^{6} r - 4 \hat{h}^{5} a k m^{5} r^{2} - 4 \hat{h}^{5} e h m^{4} r^{2} + 16 \hat{h}^{5} h m^{5} r^{2} + \hat{h}^{4} a^{2} e^{2} k^{2} m^{4} - 2 \hat{h}^{4} a^{2} e k^{2} m^{5} - 2 \hat{h}^{4} a^{2} e k^{2} m^{4} r + \hat{h}^{4} a^{2} k^{2} m^{6} + 2 \hat{h}^{4} a^{2} k^{2} m^{5} r + \hat{h}^{4} a^{2} k^{2} m^{4} r^{2} + 2 \hat{h}^{4} a e^{2} h k m^{3} r + 14 \hat{h}^{4} a e h k m^{4} r + 6 \hat{h}^{4} a e h k m^{3} r^{2} - 16 \hat{h}^{4} a h k m^{5} r - 16 \hat{h}^{4} a h k m^{4} r^{2} + \hat{h}^{4} e^{2} h^{2} m^{2} r^{2} - 8 \hat{h}^{4} e h^{2} m^{3} r^{2} + 24 \hat{h}^{4} h^{2} m^{4} r^{2} + 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{3} - 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{2} r - 6 \hat{h}^{3} a^{2} e h k^{2} m^{4} - 8 \hat{h}^{3} a^{2} e h k^{2} m^{3} r - 2 \hat{h}^{3} a^{2} e h k^{2} m^{2} r^{2} + 4 \hat{h}^{3} a^{2} h k^{2} m^{5} + 8 \hat{h}^{3} a^{2} h k^{2} m^{4} r + 4 \hat{h}^{3} a^{2} h k^{2} m^{3} r^{2} + 2 \hat{h}^{3} a e^{2} h^{2} k m^{2} r - 2 \hat{h}^{3} a e^{2} h^{2} k m r^{2} + 16 \hat{h}^{3} a e h^{2} k m^{3} r + 12 \hat{h}^{3} a e h^{2} k m^{2} r^{2} - 24 \hat{h}^{3} a h^{2} k m^{4} r - 24 \hat{h}^{3} a h^{2} k m^{3} r^{2} - 4 \hat{h}^{3} e h^{3} m^{2} r^{2} + 16 \hat{h}^{3} h^{3} m^{3} r^{2} + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m^{2} - 2 \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m r + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} r^{2} - 6 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{3} - 10 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{2} r - 4 \hat{h}^{2} a^{2} e h^{2} k^{2} m r^{2} + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{4} + 12 \hat{h}^{2} a^{2} h^{2} k^{2} m^{3} r + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{2} r^{2} + 6 \hat{h}^{2} a e h^{3} k m^{2} r + 6 \hat{h}^{2} a e h^{3} k m r^{2} - 16 \hat{h}^{2} a h^{3} k m^{3} r - 16 \hat{h}^{2} a h^{3} k m^{2} r^{2} + 4 \hat{h}^{2} h^{4} m^{2} r^{2} - 2 \hat{h} a^{2} e h^{3} k^{2} m^{2} - 4 \hat{h} a^{2} e h^{3} k^{2} m r - 2 \hat{h} a^{2} e h^{3} k^{2} r^{2} + 4 \hat{h} a^{2} h^{3} k^{2} m^{3} + 8 \hat{h} a^{2} h^{3} k^{2} m^{2} r + 4 \hat{h} a^{2} h^{3} k^{2} m r^{2} - 4 \hat{h} a h^{4} k m^{2} r - 4 \hat{h} a h^{4} k m r^{2} + a^{2} h^{4} k^{2} m^{2} + 2 a^{2} h^{4} k^{2} m r + a^{2} h^{4} k^{2} r^{2}} - \frac{1}{2 a k \left(\hat{h} m + h\right)^{2}} \left(2 \hat{h}^{3} m^{3} r - \hat{h}^{2} a e k m^{2} + \hat{h}^{2} a k m^{3} - \hat{h}^{2} a k m^{2} r - \hat{h}^{2} e h m r + 4 \hat{h}^{2} h m^{2} r - \hat{h} a e h k m + \hat{h} a e h k r + 2 \hat{h} a h k m^{2} - 2 \hat{h} a h k m r + 2 \hat{h} h^{2} m r + a h^{2} k m - a h^{2} k r\right), \quad 1, \quad \left [ \left[\begin{matrix}- \frac{2 \hat{h} a d k m \left(\hat{h}^{2} m^{2} + 2 \hat{h} h m + h^{2}\right)}{\left(\hat{h} m + h\right) \left(2 \hat{h}^{3} m^{3} r + \hat{h}^{2} a e k m^{2} - \hat{h}^{2} a k m^{3} - \hat{h}^{2} a k m^{2} r - \hat{h}^{2} e h m r + 4 \hat{h}^{2} h m^{2} r + \hat{h} a e h k m + \hat{h} a e h k r - 2 \hat{h} a h k m^{2} - 2 \hat{h} a h k m r + 2 \hat{h} h^{2} m r - a h^{2} k m - a h^{2} k r + \sqrt{4 \hat{h}^{6} m^{6} r^{2} + 4 \hat{h}^{5} a e k m^{5} r - 4 \hat{h}^{5} a k m^{6} r - 4 \hat{h}^{5} a k m^{5} r^{2} - 4 \hat{h}^{5} e h m^{4} r^{2} + 16 \hat{h}^{5} h m^{5} r^{2} + \hat{h}^{4} a^{2} e^{2} k^{2} m^{4} - 2 \hat{h}^{4} a^{2} e k^{2} m^{5} - 2 \hat{h}^{4} a^{2} e k^{2} m^{4} r + \hat{h}^{4} a^{2} k^{2} m^{6} + 2 \hat{h}^{4} a^{2} k^{2} m^{5} r + \hat{h}^{4} a^{2} k^{2} m^{4} r^{2} + 2 \hat{h}^{4} a e^{2} h k m^{3} r + 14 \hat{h}^{4} a e h k m^{4} r + 6 \hat{h}^{4} a e h k m^{3} r^{2} - 16 \hat{h}^{4} a h k m^{5} r - 16 \hat{h}^{4} a h k m^{4} r^{2} + \hat{h}^{4} e^{2} h^{2} m^{2} r^{2} - 8 \hat{h}^{4} e h^{2} m^{3} r^{2} + 24 \hat{h}^{4} h^{2} m^{4} r^{2} + 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{3} - 2 \hat{h}^{3} a^{2} e^{2} h k^{2} m^{2} r - 6 \hat{h}^{3} a^{2} e h k^{2} m^{4} - 8 \hat{h}^{3} a^{2} e h k^{2} m^{3} r - 2 \hat{h}^{3} a^{2} e h k^{2} m^{2} r^{2} + 4 \hat{h}^{3} a^{2} h k^{2} m^{5} + 8 \hat{h}^{3} a^{2} h k^{2} m^{4} r + 4 \hat{h}^{3} a^{2} h k^{2} m^{3} r^{2} + 2 \hat{h}^{3} a e^{2} h^{2} k m^{2} r - 2 \hat{h}^{3} a e^{2} h^{2} k m r^{2} + 16 \hat{h}^{3} a e h^{2} k m^{3} r + 12 \hat{h}^{3} a e h^{2} k m^{2} r^{2} - 24 \hat{h}^{3} a h^{2} k m^{4} r - 24 \hat{h}^{3} a h^{2} k m^{3} r^{2} - 4 \hat{h}^{3} e h^{3} m^{2} r^{2} + 16 \hat{h}^{3} h^{3} m^{3} r^{2} + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m^{2} - 2 \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} m r + \hat{h}^{2} a^{2} e^{2} h^{2} k^{2} r^{2} - 6 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{3} - 10 \hat{h}^{2} a^{2} e h^{2} k^{2} m^{2} r - 4 \hat{h}^{2} a^{2} e h^{2} k^{2} m r^{2} + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{4} + 12 \hat{h}^{2} a^{2} h^{2} k^{2} m^{3} r + 6 \hat{h}^{2} a^{2} h^{2} k^{2} m^{2} r^{2} + 6 \hat{h}^{2} a e h^{3} k m^{2} r + 6 \hat{h}^{2} a e h^{3} k m r^{2} - 16 \hat{h}^{2} a h^{3} k m^{3} r - 16 \hat{h}^{2} a h^{3} k m^{2} r^{2} + 4 \hat{h}^{2} h^{4} m^{2} r^{2} - 2 \hat{h} a^{2} e h^{3} k^{2} m^{2} - 4 \hat{h} a^{2} e h^{3} k^{2} m r - 2 \hat{h} a^{2} e h^{3} k^{2} r^{2} + 4 \hat{h} a^{2} h^{3} k^{2} m^{3} + 8 \hat{h} a^{2} h^{3} k^{2} m^{2} r + 4 \hat{h} a^{2} h^{3} k^{2} m r^{2} - 4 \hat{h} a h^{4} k m^{2} r - 4 \hat{h} a h^{4} k m r^{2} + a^{2} h^{4} k^{2} m^{2} + 2 a^{2} h^{4} k^{2} m r + a^{2} h^{4} k^{2} r^{2}}\right)}\\1\end{matrix}\right]\right ]\right )\right ]$$



In [ ]: