Human immune response to infectious disease

Unos modelos dinámicos.



In [1]:

    
# Para hacer experimentos numéricos importamos numpy
import numpy as np

# y biblioteca para plotear
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

# cómputo simbólico con sympy
from sympy import *
# init_printing(use_latex='matplotlib')  # en emacs
init_printing()

Linfocitos sin virus



In [2]:

    
Lamda1, Lamda2, mu1, mu2, a1, a2, b1, b2, E1, E2 = symbols('Lamda1 Lamda2 mu1 mu2 a1 a2 b1 b2 E1 E2')



In [3]:

    
dE1 = Lamda1 - mu1 * E1 + (a1 * E1 * E2)/(1 + b1 * E1 * E2)
dE1









    Out[3]:





$$\frac{E_{1} E_{2} a_{1}}{E_{1} E_{2} b_{1} + 1} - E_{1} \mu_{1} + \Lambda_{1}$$



In [4]:

    
dE2 = Lamda2 - mu2 * E2 + (a2 * E1 * E2)/(1 + b2 * E1 * E2)
dE2









    Out[4]:





$$\frac{E_{1} E_{2} a_{2}}{E_{1} E_{2} b_{2} + 1} - E_{2} \mu_{2} + \Lambda_{2}$$

Búsqueda de Equilibrios

Abuso del cómputo simbólico

Los objetos creados dE1 y dE2 contienen operaciones representadas simbólicamente que conjugan todas las variables del sistema. Imprimirlas reconstruye las ecuaciones en su notación matemática, a partir de las expresiones simbólicas de Python. Para lograrlo hay que declararle a Python cuáles símbolos componen las expresiones matemáticas. Esto permite hacer operaciones simbólicas sobre las expresiones, como resolver, derivar, sustituir.

Todas las magnitudes del modelo se representan en las ecuaciones por símbolos, como $\Lambda_{1}, \mu_{2}, a{1}$, etc.

Por otro lado, las expresiones simbólicas dE1 y dE2 se pueden convertir en funciones ejecutables, pero estos símbolos no son argumentos de las funciones dE1 y dE2, se asume que serán constantes y que las únicas variables con cambios serán E1 y E2.

En el intento que sigue de encontrar algebráicamente los equilibrios del sistema, he incluido como símbolos todas esas constantes del sistema. Mi esperanza era usar operaciones de SymPy hasta tener una expresión simbólica que volver ejecutable, y sólo entonces alimentarle valores numéricos al modelo computacional.

Es que es una mala práctica de programación incluir en el código "números mágicos". O sea: números en medio de expresiones matemáticas, sin explicación, e.g.:

a = 3.14159 * r**2

Es mejor estilo escribir

pi = 3.141592
a = pi * r**2

Pero cada símbolo es una variable para SymPy, y el intento de encontrar soluciones al sistema de ecuaciones con tantas variables genera expresiones ilegibles.



In [5]:

    
# resolver dE2 para E1, viceversa
e1 = solve(dE2, E1)[0]
e2 = solve(dE1, E2)[0]
e1, e2









    Out[5]:





$$\left ( \frac{E_{2} \mu_{2} - \Lambda_{2}}{E_{2} \left(- E_{2} b_{2} \mu_{2} + \Lambda_{2} b_{2} + a_{2}\right)}, \quad \frac{E_{1} \mu_{1} - \Lambda_{1}}{E_{1} \left(- E_{1} b_{1} \mu_{1} + \Lambda_{1} b_{1} + a_{1}\right)}\right )$$



In [6]:

    
# sustituir una de las expresiones de la celda anterior en una de las expresiones originales, resolver
static_e1 = solve(e1.subs(E2, e2), E1)
static_e1









    Out[6]:





$$\left [ \frac{1}{2 \Lambda_{2} b_{1} \mu_{1}} \left(\Lambda_{1} \Lambda_{2} b_{1} + \Lambda_{2} a_{1} - \mu_{1} \mu_{2} - \sqrt{\Lambda_{1}^{2} \Lambda_{2}^{2} b_{1}^{2} + 2 \Lambda_{1} \Lambda_{2}^{2} a_{1} b_{1} + 2 \Lambda_{1} \Lambda_{2} b_{1} \mu_{1} \mu_{2} + \Lambda_{2}^{2} a_{1}^{2} - 2 \Lambda_{2} a_{1} \mu_{1} \mu_{2} + \mu_{1}^{2} \mu_{2}^{2}}\right), \quad \frac{1}{2 \Lambda_{2} b_{1} \mu_{1}} \left(\Lambda_{1} \Lambda_{2} b_{1} + \Lambda_{2} a_{1} - \mu_{1} \mu_{2} + \sqrt{\Lambda_{1}^{2} \Lambda_{2}^{2} b_{1}^{2} + 2 \Lambda_{1} \Lambda_{2}^{2} a_{1} b_{1} + 2 \Lambda_{1} \Lambda_{2} b_{1} \mu_{1} \mu_{2} + \Lambda_{2}^{2} a_{1}^{2} - 2 \Lambda_{2} a_{1} \mu_{1} \mu_{2} + \mu_{1}^{2} \mu_{2}^{2}}\right)\right ]$$

Esto es una barbaridad, un error. A continuación cómo lo resolví.

Sustituir magnitudes al crear las expresiones simbólicas

Si declaro antes, ya no son números mágicos. Al crear las expresiones se interpolarán los valores numéricos así que las expresiones simbólicas serán simples: SymPy hace la aritmética necesaria para simplificar, lo que le era imposible con tanto símbolo.



In [7]:

    
# borro las variables previas, para quitarles lo simbólico
#del(Lamda1, Lamda2, mu1, mu2, a1, a2, b1, b2, E1, E2, dE1, dE2)

# variables simples de Python
Lamda1 = 1
Lamda2 = 1
mu1 = 1.25
mu2 = 1.25
a1 = 0.252
a2 = 0.252
b1 = 0.008
b2 = 0.008

# esta vez sólo dos símbolos
E1, E2 = symbols('E1 E2')



In [8]:

    
dE1 = Lamda1 - mu1 * E1 + (a1 * E1 * E2)/(1 + b1 * E1 * E2)
dE1









    Out[8]:





$$\frac{0.252 E_{1} E_{2}}{0.008 E_{1} E_{2} + 1} - 1.25 E_{1} + 1$$



In [9]:

    
dE2 = Lamda2 - mu2 * E2 + (a2 * E1 * E2)/(1 + b2 * E1 * E2)
dE2









    Out[9]:





$$\frac{0.252 E_{1} E_{2}}{0.008 E_{1} E_{2} + 1} - 1.25 E_{2} + 1$$

¡Mucho mejor! Ahora sí:

Búsqueda de los equilibrios



In [10]:

    
e2 = solve(dE1, E2)[0]
e2









    Out[10]:





$$\frac{25.0 \left(- 5.0 E_{1} + 4.0\right)}{E_{1} \left(E_{1} - 26.0\right)}$$



In [11]:

    
solve(dE2.subs({E2: e2}), E1)









    Out[11]:





$$\left [ 1.0, \quad 5.0, \quad 20.0\right ]$$

Jacobiana



In [12]:

    
J = symbols("J")

J = Matrix([[diff(dE1, E1), diff(dE1, E2)], 
            [diff(dE2, E1), diff(dE2, E2)]])
J









    Out[12]:





$$\left[\begin{matrix}- \frac{0.002016 E_{1} E_{2}^{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{2}}{0.008 E_{1} E_{2} + 1} - 1.25 & - \frac{0.002016 E_{1}^{2} E_{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{1}}{0.008 E_{1} E_{2} + 1}\\- \frac{0.002016 E_{1} E_{2}^{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{2}}{0.008 E_{1} E_{2} + 1} & - \frac{0.002016 E_{1}^{2} E_{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{1}}{0.008 E_{1} E_{2} + 1} - 1.25\end{matrix}\right]$$

Evaluada en los puntos de equilibrio

(1,1) un sumidero



In [13]:

    
Je1 = J.subs({E1: 1, E2:1})
Je1









    Out[13]:





$$\left[\begin{matrix}-1.00198412698413 & 0.248015873015873\\0.248015873015873 & -1.00198412698413\end{matrix}\right]$$



In [14]:

    
Je1.det(), Je1.trace()









    Out[14]:





$$\left ( 0.942460317460317, \quad -2.00396825396825\right )$$



In [15]:

    
Je1.eigenvects()









    Out[15]:





$$\left [ \left ( -1.25, \quad 1, \quad \left [ \left[\begin{matrix}-1.0\\1.0\end{matrix}\right]\right ]\right ), \quad \left ( -0.753968253968257, \quad 1, \quad \left [ \left[\begin{matrix}1.0\\1.0\end{matrix}\right]\right ]\right )\right ]$$

(5, 5) un punto de ensilladura



In [16]:

    
Je2 = J.subs({E1: 5, E2:5})
Je2









    Out[16]:





$$\left[\begin{matrix}-0.375 & 0.875\\0.875 & -0.375\end{matrix}\right]$$



In [17]:

    
Je2.det(), Je2.trace()









    Out[17]:





$$\left ( -0.625, \quad -0.75\right )$$



In [18]:

    
Je2.eigenvects(), Je2.eigenvals()









    Out[18]:





$$\left ( \left [ \left ( -1.25, \quad 1, \quad \left [ \left[\begin{matrix}-1.0\\1.0\end{matrix}\right]\right ]\right ), \quad \left ( 0.5, \quad 1, \quad \left [ \left[\begin{matrix}1.0\\1.0\end{matrix}\right]\right ]\right )\right ], \quad \left \{ - \frac{5}{4} : 1, \quad \frac{1}{2} : 1\right \}\right )$$

(20, 20) otro sumidero



In [19]:

    
Je3 = J.subs({E1: 20, E2: 20})
Je3









    Out[19]:





$$\left[\begin{matrix}-0.964285714285714 & 0.285714285714286\\0.285714285714286 & -0.964285714285714\end{matrix}\right]$$



In [20]:

    
Je3.det(), Je3.trace()









    Out[20]:





$$\left ( 0.848214285714286, \quad -1.92857142857143\right )$$



In [21]:

    
Je3.eigenvects()









    Out[21]:





$$\left [ \left ( -1.25, \quad 1, \quad \left [ \left[\begin{matrix}-1.0\\1.0\end{matrix}\right]\right ]\right ), \quad \left ( -0.678571428571429, \quad 1, \quad \left [ \left[\begin{matrix}1.0\\1.0\end{matrix}\right]\right ]\right )\right ]$$

Campos de vectores



In [22]:

    
del(Lamda1, Lamda2, mu1, mu2, a1, a2, b1, b2, E1, E2, dE1, dE2)

def dE1(E1, E2):
    return Lamda1 - mu1 * E1 + (a1 * E1 * E2)/(1 + b1*E1*E2)

def dE2(E1, E2):
    return Lamda2 - mu2 * E2 + (a2 * E1 * E2)/(1 + b2*E1*E2)



In [23]:

    
Lamda1 = 1
Lamda2 = 1
mu1 = 1.25
mu2 = 1.25
a1 = 0.252
a2 = 0.252
b1 = 0.008
b2 = 0.008

plt.figure(figsize=(16, 7))
s = 1
for eq in [(1,1), (5,5), (20,20)]:
    ax = plt.subplot(1, 3, s)
    s += 1
    ax.set_title("(%s, %s)" % eq, fontsize=13)
    
    i, j = np.meshgrid(np.linspace(eq[0] - 2,
                                   eq[0] + 2, 15),
                       np.linspace(eq[1] - 2,
                                   eq[1] + 2, 15))

    u = dE1(i, j)
    v = dE2(i, j)

    plt.quiver(i, j, u, v)

Trayectorias



In [24]:

    
def step(x, y, dt, f, g):
    return (x + dt * f(x, y),
             y + dt * g(x, y))

def trayectoria(x0, y0, f, g, dt=0.01, steps=100):
    x = x0
    y = y0
    t = list()
    for n in range(steps):
        t.append((x, y))
        x, y = step(x, y, dt, f, g)
    return t



In [25]:

    
def linea(E1_0, E2_0, dt=0.1, steps=300, color='teal'):
    E1, E2 = zip(*[coord for coord in
                   trayectoria(E1_0, E2_0, dE1, dE2, dt, steps)])
    plt.plot(E1, E2, color=color)
    plt.arrow(E1[-2], E2[-2],E1[-1] - E1[-2], E2[-1] - E2[-2],
               color=color, head_width=0.3, head_length=0.3)


i, j = np.meshgrid(np.linspace(-3, 24, 25),
                   np.linspace(-3, 24, 25))

u = dE1(i, j)
v = dE2(i, j)

fig = plt.figure(figsize=(9,7))
plt.quiver(i, j, u, v)

linea(0,0, dt=0.01)
linea(0,14, dt=0.01)
linea(10,-1, dt=0.01, steps=120)


for N in np.linspace(-0.5, 0.5, 32):
    linea(24, N, dt=0.01, color='orange')
for N in np.linspace(-0.5, 0.5, 32):
    linea(N, 24, dt=0.01, color='violet')


linea(22, 20, dt=0.01, color='purple')
linea(20, 22, dt=0.01, color='purple')
linea(10, 20, dt=0.01, color='purple')
linea(20, 6, dt=0.01, color='purple')

plt.plot(1, 1, 'teal', marker='o', markersize=10, alpha=0.4)
plt.plot(5, 5, 'orange', marker='o', markersize=10, alpha=0.4)
plt.plot(20, 20, 'purple', marker='o', markersize=10, alpha=0.4)

plt.show()

Tamaños de población al paso del tiempo



In [26]:

    
fig = plt.figure(figsize=(10, 6))

def pop(E1_0, E2_0, steps = 1500):
    e1t = [E1_0, ] + np.zeros(steps)
    e2t = [E2_0, ] + np.zeros(steps)

    dt = 0.01
    for t in range(1, steps):
        e1t[t] = e1t[t-1] + dt * dE1(e1t[t-1], e2t[t-1])
        e2t[t] = e2t[t-1] + dt * dE2(e1t[t-1], e2t[t-1])

    plt.plot(range(steps), e1t, color='teal')
    plt.plot(range(steps), e2t, color='orange')

pop(2, 3)
pop(4,6)
pop(18,22)

En presencia del Virus



In [34]:

    
# del(a1, a2, b1, b2, dE1, dE2)

Lamda1 = 1
Lamda2 = 1
mu1 = 1.25
mu2 = 1.25
a1 = 0.252
a2 = 0.252
b1 = 0.008
b2 = 0.008
r = 0.07
k = 0.01
K = 0.05

V, E1, E2 = symbols('V E1 E2')

dE1 = Lamda1 - mu1 * E1 + (a1 * E1 * E2)/(1 + b1 * E1 * E2) + K*V*E1
dE2 = Lamda2 - mu2 * E2 + (a2 * E1 * E2)/(1 + b2 * E1 * E2)
dV = r*V - k*V*E1



In [37]:

    
J = symbols("J")

J = Matrix([[diff(dE1, E1), diff(dE1, E2), diff(dE1, V)], 
            [diff(dE2, E1), diff(dE2, E2), diff(dE2, V)],
            [diff(dV, E1), diff(dV, E2), diff(dV, V)]])
J









    Out[37]:





$$\left[\begin{matrix}- \frac{0.002016 E_{1} E_{2}^{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{2}}{0.008 E_{1} E_{2} + 1} + 0.05 V - 1.25 & - \frac{0.002016 E_{1}^{2} E_{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{1}}{0.008 E_{1} E_{2} + 1} & 0.05 E_{1}\\- \frac{0.002016 E_{1} E_{2}^{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{2}}{0.008 E_{1} E_{2} + 1} & - \frac{0.002016 E_{1}^{2} E_{2}}{\left(0.008 E_{1} E_{2} + 1\right)^{2}} + \frac{0.252 E_{1}}{0.008 E_{1} E_{2} + 1} - 1.25 & 0\\- 0.01 V & 0 & - 0.01 E_{1} + 0.07\end{matrix}\right]$$



In [38]:

    
Je1 = J.subs({E1: 1, E2:1, V:0})
Je1









    Out[38]:





$$\left[\begin{matrix}-1.00198412698413 & 0.248015873015873 & 0.05\\0.248015873015873 & -1.00198412698413 & 0\\0 & 0 & 0.06\end{matrix}\right]$$



In [39]:

    
Je1 = J.subs({E1: 5, E2:5, V:0})
Je1









    Out[39]:





$$\left[\begin{matrix}-0.375 & 0.875 & 0.25\\0.875 & -0.375 & 0\\0 & 0 & 0.02\end{matrix}\right]$$



In [40]:

    
Je1 = J.subs({E1: 20, E2:20, V:0})
Je1









    Out[40]:





$$\left[\begin{matrix}-0.964285714285714 & 0.285714285714286 & 1.0\\0.285714285714286 & -0.964285714285714 & 0\\0 & 0 & -0.13\end{matrix}\right]$$



In [41]:

    
def dE1(E1, E2, V):
    return Lamda1 - mu1 * E1 + (a1 * E1 * E2)/(1 + b1*E1*E2) + K*V*E1

def dE2(E1, E2):
    return Lamda2 - mu2 * E2 + (a2 * E1 * E2)/(1 + b2*E1*E2)


def dV(V, E1):
    return r*V - k*V*E1



In [68]:

    
from mpl_toolkits.mplot3d import Axes3D

Lamda1 = 1
Lamda2 = 1
mu1 = 1.25
mu2 = 1.25
a1 = 0.252
a2 = 0.252
b1 = 0.008
b2 = 0.008
r = 0.07
k = 0.01
K = 0.05

def pops(e1_0, e2_0, v_0, steps=10000, dt=0.001):

    e1t = [e1_0, ] + np.zeros(steps)
    e2t = [e2_0, ] + np.zeros(steps)
    vt = [v_0, ] + np.zeros(steps)

    for t in range(1, steps):
        e1t[t] = e1t[t-1] + dt * dE1(e1t[t-1], e2t[t-1], vt[t-1])
        e2t[t] = e2t[t-1] + dt * dE2(e1t[t-1], e2t[t-1])
        vt[t] = vt[t-1] + dt * dV(e1t[t-1], vt[t-1])
    return e1t, e2t, vt


fig = plt.figure(figsize=(14, 13))
ax = fig.gca(projection='3d')
for i in np.linspace(0, 22, 10):
    for j in np.linspace(0, 22, 10):
        e1t, e2t, vt = pops(i, j, 0.1)
        ax.plot(vt, e1t, e2t)

ax.set_xlabel("V")
ax.set_ylabel("$E_{1}$")
ax.set_zlabel("$E_{2}$")
plt.show()



In [69]:

    
fig = plt.figure(figsize=(10,6))
e1t, e2t, vt = pops(20, 20, 30)
plt.plot(e1t, color='teal')
plt.plot(e2t, color='orange')
plt.plot(vt, color='green')









    Out[69]:





[<matplotlib.lines.Line2D at 0x7fa0a1a93c88>]



In [70]:

    
fig = plt.figure(figsize=(10,6))
e1t, e2t, vt = pops(20, 10, 30)
plt.plot(e1t, color='teal')
plt.plot(e2t, color='orange')
plt.plot(vt, color='green')









    Out[70]:





[<matplotlib.lines.Line2D at 0x7fa09f7b0748>]



In [71]:

    
fig = plt.figure(figsize=(10,6))
e1t, e2t, vt = pops(2, 1, 30)
plt.plot(e1t, color='teal')
plt.plot(e2t, color='orange')
plt.plot(vt, color='green')









    Out[71]:





[<matplotlib.lines.Line2D at 0x7fa09f7e0160>]



In [74]:

    
fig = plt.figure(figsize=(10,6))
e1t, e2t, vt = pops(1, 1, 0.1, steps=60000)
plt.plot(e1t, color='teal')
plt.plot(e2t, color='orange')
plt.plot(vt, color='green')









    Out[74]:





[<matplotlib.lines.Line2D at 0x7fa09f7cdf60>]