Infectious Pulse

https://github.com/alvason/infectious-pulse/

Many-strain SIR evolution --- vaccine events


In [6]:
'''
author: Alvason Zhenhua Li
date:   03/23/2015
'''
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
import time

import alva_machinery_event as alva

AlvaFontSize = 23
AlvaFigSize = (12, 7)
numberingFig = 0

numberingFig = numberingFig + 1
plt.figure(numberingFig, figsize=(12, 6))
plt.axis('off')
plt.title(r'$ n-strain \ SIRV \ equations \ (pulsed \ vaccination) $',fontsize = AlvaFontSize)
plt.text(0, 6.0/8,r'$ \frac{\partial S_i(t)}{\partial t} = \
        + \mu N - \mu S_i(t) - \beta S_i(t)\sum_{j = 1}^{n} (1 - \frac{|j - i|}{r + |j - i|})I_{i}(t) $'
         , fontsize = 1.2*AlvaFontSize)
plt.text(0, 4.0/8, r'$ \frac{\partial I_i(t)}{\partial t} = \
         + \beta S_i(t)I_i(t) - \gamma I_i(t) - \mu I_i(t) \
         + m \frac{I_{i - 1}(t) - 2I_i(t) + I_{i + 1}(t)}{(\Delta i)^2} $'
         , fontsize = 1.2*AlvaFontSize)
plt.text(0, 2.0/8, r'$ \frac{\partial V_i(t)}{\partial t} = \
         + \phi S_i(t) \xi (i - k, t - \tau, l) - \gamma V_i(t) - \mu V_i(t) $'
         , fontsize = 1.2*AlvaFontSize)
plt.text(0, 0.0/8,r'$ \frac{\partial R_i(t)}{\partial t} = \
         +\gamma I_i(t) - \mu R_i(t) - \beta S_i(t)I_i(t) \
         + \beta S_i(t)\sum_{j = 1}^{n} (1 - \frac{|j - i|}{r + |j - i|})I_{i}(t) \
         + \phi S_i(t) \xi (i - k, t - \tau, l)\sum_{j = 1}^{n} (1 - \frac{|j - i|}{r + |j - i|}) $'
         , fontsize = 1.2*AlvaFontSize)
plt.show()

# define many-strain S-I-R equation
def dSdt_array(SIRxt = [], *args):
    # naming
    S = SIRxt[0]
    I = SIRxt[1]
    R = SIRxt[2]
    x_totalPoint = SIRxt.shape[1]
    # there are n dRdt
    dS_dt_array = np.zeros(x_totalPoint)
    # each dSdt with the same equation form
    dS_dt_array[:] = - infecRate*S[:]*crossI_neighborSum_X(I, cross_radius, gX)[:] \
                     + inOutRate*totalSIR - inOutRate*S[:]
    return(dS_dt_array)

def dIdt_array(SIRxt = [], *args):
    global mutatRate
    mutatRate = alva.eventName
    # naming
    S = SIRxt[0]
    I = SIRxt[1]
    R = SIRxt[2]
    x_totalPoint = SIRxt.shape[1]
    # there are n dRdt
    dI_dt_array = np.zeros(x_totalPoint)
    # each dIdt with the same equation form
    Icopy = np.copy(I)
    centerX = Icopy[:]
    leftX = np.roll(Icopy[:], 1)
    rightX = np.roll(Icopy[:], -1)
    leftX[0] =centerX[0]
    rightX[-1] = centerX[-1]
    dI_dt_array[:] = + infecRate*S[:]*I[:] \
                     - recovRate*I[:] - inOutRate*I[:] \
                     + mutatRate*(leftX[:] - 2*centerX[:] + rightX[:])/(dx**2)                                                                                             
    return(dI_dt_array)

def dRdt_array(SIRxt = [], *args):
    # naming
    S = SIRxt[0]
    I = SIRxt[1]
    R = SIRxt[2]
    x_totalPoint = SIRxt.shape[1]
    # there are n dRdt
    dR_dt_array = np.zeros(x_totalPoint)
    # each dIdt with the same equation form
    dR_dt_array[:] = + recovRate*I[:] - inOutRate*R[:] \
                     - infecRate*S[:]*I[:] \
                     + infecRate*S[:]*crossI_neighborSum_X(I, cross_radius, gX)[:]
    return(dR_dt_array)

# inverted-monod equation
def monodInvert(half_radius, i):
    if half_radius == 0:
        gOut = i*0
        # numpy.reshape will not change the structure of i, 
        # so that the first element of i(unkonwn-size-array) can be setted by array_to_list[0] 
        array_to_list = np.reshape(i,[i.size,1]) 
        array_to_list[0] = 1 
    else: gOut = 1 - np.absolute(i)/(half_radius + np.absolute(i))
    return (gOut)

# cross immunity
def crossI_neighborSum_X(gI, half_radius, gX):
    total_neighbor_X = gX.shape[0]
    I_neighborSum = np.zeros(total_neighbor_X)
    # all I[xn] with neighbor-sum 
    ratioM = np.zeros([total_neighbor_X, total_neighbor_X])
    gXX = np.tile(gX, [total_neighbor_X, 1])
    gII = np.tile(gI, [total_neighbor_X, 1])
    ratioM[:, :] = monodInvert(half_radius, gXX[:, :] - gXX[:, :].T)
    I_neighborSum[:] = np.sum(ratioM[:, :] * gII[:, :].T, axis = 0)
    if half_radius == 0:
        I_neighborSum = np.copy(gI)
    return (I_neighborSum)



In [2]:
# setting parameter
timeUnit = 'year'
if timeUnit == 'day':
    day = 1
    year = 365
elif timeUnit == 'year':
    year = 1
    day = float(1)/365 
    
totalSIR = float(1) # total population
reprodNum = 1.8 # basic reproductive number R0: one infected person will transmit to 1.8 person 
recovRate = float(1)/(4*day) # 4 days per period ==> rate/year = 365/4
inOutRate = float(1)/(30*year) # birth rate per year
infecRate = reprodNum*(recovRate + inOutRate)/totalSIR # per year, per person, per total-population
# mutatRate = float(1)/(10**17) # mutation rate
cross_radius = float(0) # radius of cross-immunity (the distance of half-of-value in the Monod equation) 

# time boundary and griding condition
minT = float(0)*year
maxT = float(5)*year
totalPoint_T = int(1*10**3 + 1)
spacingT = np.linspace(minT, maxT, num = totalPoint_T, retstep = True)
gT = spacingT[0]
dt = spacingT[1]

# space boundary and griding condition
minX = float(0)
maxX = float(15)
totalPoint_X = int(maxX*1 + 1)
gridingX = np.linspace(minX, maxX, num = totalPoint_X, retstep = True)
gX = gridingX[0]
dx = gridingX[1]

gS_array = np.zeros([totalPoint_X, totalPoint_T])
gI_array = np.zeros([totalPoint_X, totalPoint_T])
gR_array = np.zeros([totalPoint_X, totalPoint_T])

# initial output condition (only one virus in equilibrium condition)

gI_array[0, 0] = float(1)/10**6  # only one virus exists
gR_array[0, 0] = float(0)
gS_array[:, 0] = totalSIR - gI_array[:, 0] - gR_array[:, 0]

event_tn_In = np.array([[0, float(1)/(10**17)], [2, float(1)/(10**7)]]) # cross_radius 

# Runge Kutta numerical solution
time_start = time.time()
pde_array = np.array([dSdt_array, dIdt_array, dRdt_array])
initial_Out = np.array([gS_array, gI_array, gR_array])
gOut_array = alva.AlvaRungeKutta4XT(pde_array, initial_Out, minX, maxX, totalPoint_X, minT, maxT, totalPoint_T, event_tn_In)
time_stop = time.time()
print ('machinery running time =', time_stop - time_start)
# plotting
gS = gOut_array[0]  
gI = gOut_array[1]
gR = gOut_array[2]

numberingFig = numberingFig + 1
maxLevel = gI_array[0, 0]*10**5.2
plt.figure(numberingFig, figsize = AlvaFigSize)
plt.contourf(gT, gX, gI, levels = np.arange(0, maxLevel, maxLevel/100))
plt.title(r'$ Infectious \ pulse \ by \ mutation \ and \ cross-immunity $', fontsize = AlvaFontSize)
plt.xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
plt.ylabel(r'$ discrete \ space \ (strain) $', fontsize = AlvaFontSize)
plt.colorbar()
plt.text(maxT*4.0/3, maxX*5.0/6, r'$ R_0 = %f $'%(reprodNum), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*4.0/6, r'$ \gamma = %f $'%(recovRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*3.0/6, r'$ \beta = %f $'%(infecRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*2.0/6, r'$ \mu = %f $'%(inOutRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*1.0/6, r'$ m = %f, \ \ m = %f $'%(event_tn_In[0, 1]*10**14, event_tn_In[1, 1]*10**4*5), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*0.0/6, r'$ r = %f $'%(cross_radius), fontsize = AlvaFontSize)
plt.show()


('machinery running time =', 9.895304203033447)

In [3]:
# Normalization stacked graph
cutting_X = -1
cutting_T = -1
gI_N = np.copy(gI)
for xn in range(totalPoint_X):
    gI_N[xn, :] = gI_N[xn, :]/np.sum(gI_N[xn, :]*dt)
numberingFig = numberingFig + 1
plt.figure(numberingFig, figsize = (16, 4))
plt.stackplot(gT[:cutting_T], gI_N[:cutting_X, :cutting_T], alpha = 0.3)
plt.title(r'$ Normalization-stacked-graph \ of \ infectious \ pulse---Prevalence $', fontsize = AlvaFontSize)
plt.xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
plt.ylabel(r'$ I(n,t) $', fontsize = AlvaFontSize)
plt.xlim(minT, maxT - 1)
#plt.ylim(0, 4)
plt.axes().set_xticks(np.arange(minT, maxT))
plt.grid(True)
plt.show()

# Proportion stacked graph
gI_P = np.copy(gI)
for tn in range(totalPoint_T):
    gI_P[:, tn] = gI_P[:, tn]/np.sum(gI_P[:, tn])
numberingFig = numberingFig + 1
plt.figure(numberingFig, figsize = (16, 4))
plt.stackplot(gT[:cutting_T], gI_P[:cutting_X, :cutting_T], alpha = 0.3)
plt.title(r'$ Proportion-stacked-graph \ of \ infectious \ pulse---Prevalence $', fontsize = AlvaFontSize)
plt.xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
plt.ylabel(r'$ I(n,t) $', fontsize = AlvaFontSize)
plt.xlim(minT, maxT - 1)
#plt.ylim(0, 1)
plt.axes().set_xticks(np.arange(minT, maxT))
plt.grid(True)
plt.show()

# Normalization stacked graph of Incidence
gInc = np.zeros([totalPoint_X, totalPoint_T])
for xn in range(totalPoint_X):
    gInc[xn, :] = infecRate*gS[xn, :]*gI[xn, :]/np.sum(infecRate*gS[xn, :]*gI[xn, :]*dt)
numberingFig = numberingFig + 1
plt.figure(numberingFig, figsize = (16, 4))
plt.stackplot(gT[:cutting_T], gInc[:cutting_X, :cutting_T], alpha = 0.3)
plt.title(r'$ Normalization-stacked-graph \ of \ infectious \ pulse---Incidence $', fontsize = AlvaFontSize)
plt.xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
plt.ylabel(r'$ \beta S(n,t)I(n,t) $', fontsize = AlvaFontSize)
plt.xlim(minT, maxT - 1)
#plt.ylim(0, 4)
plt.axes().set_xticks(np.arange(minT, maxT))
plt.grid(True)
plt.show()



In [4]:
# plot by listing each strain 
numberingFig = numberingFig + 1
for i in range(0, totalPoint_X, 10):
    figure = plt.figure(numberingFig, figsize = (12, 3))
    plot1 = figure.add_subplot(1, 1, 1)
    plot1.plot(gT, gS[i], label = r'$ S_{%i}(t) $'%(i), color = 'blue')
    plot1.plot(gT, gR[i], label = r'$ R_{%i}(t) $'%(i), color = 'green')
    plot1.plot(gT, (gS[i] + gI[i] + gR[i]).T, label = r'$ S_{%i}(t)+I_{%i}(t)+R_{%i}(t) $'%(i, i, i), color = 'black')    
    plot1.set_xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
    plot1.set_ylabel(r'$ S \ , \ R $', fontsize = AlvaFontSize)
    plot1.set_ylim(0, totalSIR*1.1)
    plot1.legend(loc = (1.1, 0))
    
    plot2 = plot1.twinx()
    plot2.plot(gT, gI[i], label = r'$ I_{%i}(t) $'%(i), color = 'red')
    plot2.set_ylabel(r'$ I $', fontsize = AlvaFontSize, color = 'red')
    for tl in plot2.get_yticklabels(): tl.set_color('red')
    plot2.legend(loc = (1.1, 0.5))
    
    plt.grid(True)
    plt.title(r'$ SIR \ of \ strain-%i $'%(i), fontsize = AlvaFontSize)
    plt.show()



In [5]:
# 3D plotting
# define GridXX function for making 2D-grid from 1D-grid
def AlvaGridXX(gX, totalPoint_Y):
    gXX = gX;
    for n in range(totalPoint_Y - 1):
        gXX = np.vstack((gXX, gX))
    return (gXX)
# for 3D plotting
X = AlvaGridXX(gT, totalPoint_X)
Y = AlvaGridXX(gX, totalPoint_T).T 
Z = gI
numberingFig = numberingFig + 1
figure = plt.figure(numberingFig, figsize=(16, 7))

plot1 = figure.add_subplot(1, 2, 1, projection = '3d')
plot1.view_init(30, -80)
plot1.plot_wireframe(X, Y, Z, cstride = totalPoint_T, rstride = int(dx), alpha = 0.6, color = 'red')
plt.xlabel(r'$t \ (time)$', fontsize = AlvaFontSize)
plt.ylabel(r'$x \ (virus \ space)$', fontsize = AlvaFontSize)
plt.xlim(minT, maxT)
plt.ylim(minX, maxX)

plot2 = figure.add_subplot(1, 2, 2, projection = '3d')
plot2.view_init(30, 10)
plot2.plot_wireframe(X, Y, Z, cstride = totalPoint_T/40, rstride = int(maxX), alpha = 0.6)
plt.xlabel(r'$t \ (time)$', fontsize = AlvaFontSize)
plt.ylabel(r'$x \ (virus \ space)$', fontsize = AlvaFontSize)
plt.xlim(minT, maxT)
plt.ylim(minX, maxX)

figure.tight_layout()
plt.show()



In [6]:
numberingFig = numberingFig + 1;
figure = plt.figure(numberingFig, figsize = (14, 7))
plot1 = figure.add_subplot(1, 2, 2)
#colorPlot = plot1.contourf(gT, gX, gI, levels = np.arange(0, 0.12, 0.001))
colorPlot = plot1.pcolor(gT, gX, gI, vmin = 0, vmax = 0.001)
plt.title(r'$ Infectious \ pulse \ I(x,t) \ by \ mutation $', fontsize = AlvaFontSize)
plt.xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
plt.ylabel(r'$ discrete \ space \ (strain) $', fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*5.0/6, r'$ R_0 = %f $'%(reprodNum), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*4.0/6, r'$ \gamma = %f $'%(recovRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*3.0/6, r'$ \beta = %f $'%(infecRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*2.0/6, r'$ \mu = %f $'%(inOutRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*1.0/6, r'$ m = %f $'%(mutatRate*10**14), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*0.0/6, r'$ r = %i, \ \ r = %i $'%(event_tn_In[0, 1], event_tn_In[1, 1]), fontsize = AlvaFontSize)
plt.colorbar(colorPlot)
plt.xlim(minT, maxT)
plt.ylim(minX, maxX)

plot2 = figure.add_subplot(1, 2, 1)
plot2.plot(gT, (gI.T + np.arange(totalPoint_X)*0.005))
plt.title(r'$ Infectious \ pulse \ I(x,t) \ by \ mutation  $', fontsize = AlvaFontSize)
plt.xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
plt.ylabel(r'$ discrete \ space \ (strain) $', fontsize = AlvaFontSize)
plt.xlim(minT, maxT)

figure.tight_layout()
plt.show()



In [7]:
numberingFig = numberingFig + 1
maxLevel = gI_array[0, 0]*200000
plt.figure(numberingFig, figsize = AlvaFigSize)
plt.contourf(gT, gX, gI, levels = np.arange(0, maxLevel, maxLevel/100))
plt.title(r'$ Infectious \ pulse \ by \ mutation \ and \ cross-immunity $', fontsize = AlvaFontSize)
plt.xlabel(r'$time \ (%s)$'%(timeUnit), fontsize = AlvaFontSize)
plt.ylabel(r'$ discrete \ space \ (strain) $', fontsize = AlvaFontSize)
plt.colorbar()
plt.text(maxT*4.0/3, maxX*5.0/6, r'$ R_0 = %f $'%(reprodNum), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*4.0/6, r'$ \gamma = %f $'%(recovRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*3.0/6, r'$ \beta = %f $'%(infecRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*2.0/6, r'$ \mu = %f $'%(inOutRate), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*1.0/6, r'$ m = %f $'%(mutatRate*10**14), fontsize = AlvaFontSize)
plt.text(maxT*4.0/3, maxX*0.0/6, r'$ r = %f, \ \ r = %f $'%(event_tn_In[0, 1], event_tn_In[1, 1]), fontsize = AlvaFontSize)
plt.show()



In [7]: