Hard Problem

Now we turn to a harder problem that does not have a one to one mapping between all the transitions and the terms in the ODEs. We use the model in Influenza_SLIARN(), defined by

\[\begin{split}\frac{dS}{dt} &= -S \beta (I + \delta A) \\ \frac{dL}{dt} &= S \beta (I + \delta A) - \kappa L \\ \frac{dI}{dt} &= p \kappa L - \alpha I \\ \frac{dA}{dt} &= (1 - p) \kappa L - \eta A \\ \frac{dR}{dt} &= f \alpha I + \eta A \\ \frac{dN}{dt} &= -(1 - f) \alpha I.\end{split}\]

The outflow of state L, \(\kappa L\), is composed of two transitions, one to I and the other to A but the ode of L only reflects the total flow going out of the state. Same can be said for state I, where the flow \(\alpha I\) goes to both R and N. Graphically, it is a rather simple process as shown below.

We slightly change the model by introducing a new state D to convert it into a closed system. The combination of state D and N is a constant, the total population. So we can remove N and this new system consist of six transitions. We define them explicitly as ODEs and unroll them into transitions.

In [1]: from pygom import SimulateOde, Transition, TransitionType

In [2]: stateList = ['S', 'L', 'I', 'A', 'R', 'D']

In [3]: paramList = ['beta', 'p', 'kappa', 'alpha', 'f', 'delta', 'epsilon', 'N']

In [4]: odeList = [
   ...:            Transition(origin='S', equation='- beta*S/N*(I + delta*A)', transition_type=TransitionType.ODE),
   ...:            Transition(origin='L', equation='beta*S/N*(I + delta*A) - kappa*L', transition_type=TransitionType.ODE),
   ...:            Transition(origin='I', equation='p*kappa*L - alpha*I', transition_type=TransitionType.ODE),
   ...:            Transition(origin='A', equation='(1 - p)*kappa * L - epsilon*A', transition_type=TransitionType.ODE),
   ...:            Transition(origin='R', equation='f*alpha*I + epsilon*A', transition_type=TransitionType.ODE),
   ...:            Transition(origin='D', equation='(1 - f)*alpha*I', transition_type=TransitionType.ODE) ]
   ...: 

In [5]: ode = SimulateOde(stateList, paramList, ode=odeList)

In [6]: ode.get_transition_matrix()
Out[6]: 
Matrix([
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0]])

In [7]: ode2 = ode.get_unrolled_obj()

In [8]: ode2.get_transition_matrix()
Out[8]: 
Matrix([
[0, A*S*beta*delta/N + I*S*beta/N,         0,                    0,         0,                    0],
[0,                             0, L*kappa*p, -L*kappa*p + L*kappa,         0,                    0],
[0,                             0,         0,                    0, I*alpha*f, -I*alpha*f + I*alpha],
[0,                             0,         0,                    0, A*epsilon,                    0],
[0,                             0,         0,                    0,         0,                    0],
[0,                             0,         0,                    0,         0,                    0]])

In [9]: ode2.get_ode_eqn()
Out[9]: 
Matrix([
[         -A*S*beta*delta/N - I*S*beta/N],
[A*S*beta*delta/N + I*S*beta/N - L*kappa],
[                   -I*alpha + L*kappa*p],
[       -A*epsilon - L*kappa*p + L*kappa],
[                  A*epsilon + I*alpha*f],
[                   -I*alpha*f + I*alpha]])

After unrolling the odes, we have the following transition graph

In [10]: ode2.get_transition_graph()
Out[10]: <graphviz.dot.Digraph at 0x7f7568b1b710>

In [11]: plt.close()

In [12]: print(sum(ode.get_ode_eqn() - ode2.get_ode_eqn()).simplify()) # difference
0

which is exactly the same apart from slightly weird arrangement of symbols in some of the equations. The last line with a value of zero also reaffirms the result.