Estimating a parameter in an ODE at a certain time point, given other conditions

Question

Assume if I have all but one parameters in my ODE system. And I wish to infer this. Would I have to simply rearrange the equation to isolate the value? How is that done in a system where you have several equations? For example:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

#three compartments, Susceptible S, infected I, recovered R
#dS/dt, dI/dt, dR/dt
#susceptible sees birth rate coming in, deaths leaving and force of infection leaving
#infected sees FOI coming in, deaths leaving and recovery rates
#recovered sees recovery rate coming in, deaths leaving
#beta is tranmission coefficient, FOI is beta * (I/N) where N is total pop
#initially consider a model not accounting for births and deaths




# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0 = 10, 0
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0 - R0
J0 = I0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
#beta = 
gamma = 1/7
# A grid of time points (in days)
t = np.linspace(0, 100, 100+1)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I, R, J = y
    dS = ((-beta * S * I) / N)
    dI = ((beta * S * I) / N) - (gamma * I)
    dR = (gamma * I)
    dJ = ((beta * S * I) / N)
    return dS, dI, dR, dJ

# Initial conditions are S0, I0, R0
# Integrate the SIR equations over the time grid, t.
solve = odeint(deriv, (S0, I0, R0, J0), t, args=(N, beta, gamma))
S, I, R, J = solve.T

As you can see, beta I have left empty, commented out. If I have all the other values, and know that at the peak of the epidemic, 10% of the population is infected, can beta be found from all the information? What I tried was this:

sol= solve_ivp(lambda beta: deriv,
               [t], t_eval= t)
print(sol)

However this syntax does not work, I have realised. What is wrong about my approach? How can I estimate a value for beta?

Answer 1

The easiest approach here is to parameterize your code above by beta and plot the result, which is peak infections for you, as a function of beta, and then see where it crosses the treshold. Define the function:

def peak_infections_pct(beta, n_days_total = 100):

    # Total population, N.
    N = 1000
    # Initial number of infected and recovered individuals, I0 and R0.
    I0, R0 = 10, 0
    # Everyone else, S0, is susceptible to infection initially.
    S0 = N - I0 - R0
    J0 = I0
    # Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
    gamma = 1/7
    # A grid of time points (in days)
    t = np.linspace(0, n_days_total, n_days_total+1)

    # The SIR model differential equations.
    def deriv(y, t, N, beta, gamma):
        S, I, R, J = y
        dS = ((-beta * S * I) / N)
        dI = ((beta * S * I) / N) - (gamma * I)
        dR = (gamma * I)
        dJ = ((beta * S * I) / N)
        return dS, dI, dR, dJ

    # Initial conditions are S0, I0, R0
    # Integrate the SIR equations over the time grid, t.
    solve = odeint(deriv, (S0, I0, R0, J0), t, args=(N, beta, gamma))
    S, I, R, J = solve.T

    return np.max(I)/N

calculate and plot:

betas = np.linspace(0,1,101,endpoint = True)
peak_inf = [peak_infections_pct(b) for b in betas]
plt.plot(betas, peak_inf)
plt.plot(betas, 0.1*np.ones(len(betas)))

to get

so the answer is about beta ~ 0.25 To be more precise just solve for beta:

from  scipy.optimize import root
root(lambda b: peak_infections_pct(b)-0.1, x0 = 0.5).x

output:

array([0.23847079])

Note I left the time interval as an input to the function -- you may want to use different length as the epidemic may last longer that 100 days

Just to double check let's plot infections as a function of time for our beta=0.2384..:

indeed the peak is at 100 (with is 10%)