In this notebook we are going to work out the way of using the package dalgebra to the model depicted in the Flour-Beetle model. This model try to understand the evolution of a population of beetles inside the flour, were we count in a discrete way the number of larvae $L_t$, pupal $P_t$ and adult $A_t$ population with the following equations:
In this model we have two types of parameters: the $c$ family ($c_{ea}, c_{el}$ and $c_{pa}$), which appears in the exponents of the exponential functions involved; the parameter $b$ as a multiplier on how the adults put more eggs and the $\mu$ parameters ($\mu_l$ and $\mu_a$), that appear in the $(1-\mu)$ expressions throughout the linear relations between variables.
In total, we get 6 parameters. The natural question of identifiability arises in this case for different output variables.
However, we have a crucial issue in this model: the exponential functions. They appear in 2 out of 3 equations and involve 3 out of the 6 parameters. How we represent these exponential will be crucial to understand how to model this system into the dalgebra package.
import sys; sys.path.insert(0, "..") # dalgebra is here
%display latex
from dalgebra import *
R.<c_ea, c_el, c_pa, b, mu_l, mu_a> = QQ[]; R # ring for the parameters
The first idea to remove the exponentials is to substitute them by a new variable and see how this variables behaves. Let do so by: $$X_t^{(1)} = \exp\left(-c_{ea}A_t - c_{el}L_t\right),$$ $$X_t^{(2)} = \exp\left(-c_{pa}A_t\right).$$
If we try to apply the shift operator $t \mapsto t+1$ and use the equations for $A_t$ and $L_t$ we do not get anything good since some double exponentials start appearing. This is why we will (formally) consider a new operator: the derivative w.r.t. $t$. We will assume that this operators commutes with the shift we already have.
Using this derivation, we obtain the following equations:
$$X_t^{(1)'} = -(c_{ea}A_t' + c_{el}L_t') X_t^{(1)}.$$$$X_t^{(2)'} = -c_{pa}A_t'X_t^{(2)}.$$Thanks to this equations, we can eliminate the exponentials that appear in the original system but, at the same time, we will have to work on a difference-differential ring $(R, \sigma, \partial)$, which "complicates" the extension of the system we need to take to guarantee we can eliminate variables.
If we compute the necessary derivatives of our original system we obtain the following system: $$\left\{\begin{array}{ll} L_{t+1}' & = \left(b A_{t} X_t^{(1)}\right)' = \left(bA_t' - A_t(c_{ea}A_t' + c_{el}L_t')\right)X_t^{(1)} = \displaystyle\frac{bL_{t+1}A_t' - c_{ea}L_{t+1}A_tA_t' - c_{el}L_{t+1}L_t'A_t}{bA_t}\\ P_{t+1} & = (1 - \mu_l)L_t,\\ A_{t+1}' & = \left(P_t X_t^{(2)} + (1-\mu_a)A_t\right)' = \left(P_t'-c_{pa}A_t'\right)X_t^{(2)} + (1-\mu_a)A_t' = \displaystyle\frac{(P_t'- c_{pa}A_t')(A_{t+1} - (1-\mu_a)A_t)}{P_t} + (1-\mu_a)A_t' \end{array}\right.$$
The variable $P_t$ is the only variable we did not need to differentiate to remove exponentials. Moreover, its equation is so simple we can compute things like:
$$P_{t+k+1}^{(n)} = (1-\mu_l)L_{t+k}^{(n)}.$$This allows us to remove completely the variable $L_t$ in the system, since we can always represent it in terms of the $P_t$, obtaining the new system:
$$\left\{\begin{array}{ll} P_{t+2}' & = \displaystyle\frac{(1-\mu_l)(bP_{t+2}A_t' - c_{ea}P_{t+2}A_tA_t' - c_{el}P_{t+2}P_{t+1}'A_t)}{bA_t}\\ A_{t+1}' & = \displaystyle\frac{(P_t'- c_{pa}A_t')(A_{t+1} - (1-\mu_a)A_t)}{P_t} + (1-\mu_a)A_t' \end{array}\right.$$At this point we have 2 equations with 2 functions and 2 operators. We can then consider when, for example, is $A_t$ a nice output variable (i.e., whether $P_t$ can be eliminated from the system).
dalgebra¶Now it is a good moment for trying to put the system into the package dalgebra. In the current version of the code, we are not allow to have a difference-differential ring. So we will need to do it in two steps, defining more variables for the shifts of the $A_t$ and $P_t$ and then trying to get the differential structure behind.
DSR = DifferenceRing(DifferentialRing(R, lambda p : 0),R.Hom(R).one()) # we create the differential structure around these parameters
D.<A, P> = RWOPolynomialRing(DSR); D
system = RWOSystem([
(1-mu_l)*(b*P[0,2]*A[1,0] - c_ea*P[0,2]*A[0,0]*A[1,0] - c_el*P[0,2]*P[1,1]*A[0,0]) - b*A[0,0]*P[1,2],
(P[1,0]-c_pa*A[1,0])*(A[1,0]-(1-mu_a)*A[0,0]) + (1-mu_a)*A[1,0]*P[0,0] - P[0,0]*A[1,1]
], variables=[A])
system
try:
system.diff_resultant(operation = 0, alg_res="macaulay") # w.r.t. derivation
except TypeError:
print("We couldn't eliminate the As. Trying to eliminate the Ps...")
try:
system.change_variables([P]).diff_resultant(operation = 0, alg_res="macaulay") # w.r.t. derivation
except TypeError:
print("We couldn't eliminate the Ps either.")
We couldn't eliminate the As. Trying to eliminate the Ps... We couldn't eliminate the Ps either.
len(system.extend_by_operation([2,2,2], operation=0).algebraic_variables())
new_system = system.extend_by_operation([2,2,2], operation=0)
new_system
system