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

Solution 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%)