Sizing the needed resources using Python

"Flatten the curve" is one of most heard terms nowadays. Its meaning is to adopt preventive actions in order to smooth the epidemic curve avoiding to overwhelm the health system and giving enough time to prepare the system to save lives. The gentle curve is then presented against a straight line that represents the health system capacity.

In fact, the required resources are not so direct as this line. We need to consider that a patient leaves the hospital after treatment freeing resources that will be available again. Here in this article we show how to estimate the needed resources considering that, using discrete events simulation.

First part - The epidemic curve

The first step in this process is to generate the epidemic curve. A description of the mathematical model and how to generate the curve using Python is available here:

We will use this approach to generate our curve and I will show how to obtain the data to start. Suppose you have a city which population is 1.000.000 habitants and you know the number of confirmed cases of infection day by day after some weeks. Those will be our initial data. Let's consider the following:

N = 1000000
hist = pd.Series([5, 5, 5, 5, 7, 9, 14, 17, 27, 29, 31, 40, 46, 56, 66, 71])

'N' is the population and 'hist' is a Pandas Series (we will be using the Pandas library) that represents the historical number of total confirmed cases of infecteds day by day.

The first thing we could observe is that this evolution is exponential and it is possible to adjust a curve that has the following equation:

The trick here is that this exponential curve matchs the initial of the epidemic curve and we can consider the exponential growth rate as:

With this in mind we are able to know all the initial parameters to calculate the epidemic curve remembering that gamma factor for COVID-19 is considered 0,5 according to studies like [2]. Then we have:

# historical data
hist = pd.Series([5, 5, 5, 5, 7, 9, 14, 17, 27, 29, 31, 40, 46, 56, 66, 71])

# Total population N, gamma
N = 1000000
gamma = 0.5

# Exponential adjusting of the historical data
popt_exponential, pcov_exponential = \
    scipy.optimize.curve_fit(lambda t, a, b: a * np.exp(b * t),
                             hist.index, hist, p0=(5, 0.3))

# Parameters for the SIR model
I0 = popt_exponential[0]
beta = gamma + popt_exponential[1]
R0 = beta / gamma
S0 = N - I0 - R0

Now, using the epidemic model described in the previous link:

# Diferential equations of the SIR model
def deriv(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt


# Initial conditions vector
y0 = S0, I0, R0

ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T

Using the Matplotlib library and following the SciPy example (previous link) it is possible to plot our epidemic curves according to the SIR model:

Second part - The simulation

The infected I(t) is the curve to be flatten and normally is presented against a straight line that means the healthcare system capacity, see the example:

According to Brigham and Women’s Hospital in Boston, USA, 20% of infected will develop serious illness requiring hospital beds and special care. So, we will consider 20% of the infected as our red curve. Considering as example a number of 7000 beds (blue line) we see that part of the red curve is above of the blue line indicating lack of healthcare capacity, but let us see from another point of view.

Instead of taking this approach we are going to simulate the healthcare system dynamics according to the new confirmed infected cases using data from the above I(t) curve and consider the time each patient will need the hospital bed for treatment. We will use the SimPy library for the simulation. SimPy (https://simpy.readthedocs.io/en/latest/) is a process-based discrete-event simulation framework based on standard Python. With this approach is possible to simulate real life situations and is extremely useful to size required resources in many activities like logistics, manufacturing processes, airport flow, etc.

In SimPy we define processes and resources. We will define our hospital beds as limited resources and the patient's internment as our process. Let's start:

RANDOM_SEED = 100
NUM_BEDS = 7000
INTERN_TIME = random.triangular(15, 20)
SIM_TIME = 140

sys.stdout = open("./sim_dynamics.txt", "w")

'RANDOM_SEED' is a seed for our simulation, it warrants that for every run of the simulation we will get the same results, allowing us to compare different scenarios. 'NUM_BEDS' is our homecare system capacity, 'INTERN_TIME' is the time needed for the patient to be enterned using hospital bed. In this case we adopt a normal variation between 15 and 20 days. 'SIM_TIME' is the duration of the simulation in days. We also open a text file to save all prints from this point on. In this case in the same folder as the program, but you can consider a different path in your computer.

First we create a class Healthcare where we define our beds as a limited resource and also define the internment as our process which lasts the duration defined according to INTERN_TIME variable.

class Healthcare(object):
    def __init__(self, env, num_beds, intern_time):
        self.env = env
        self.bed = simpy.Resource(env, num_beds)
        self.intern_time = intern_time

    def internment(self, patient):
        yield self.env.timeout(INTERN_TIME)

Next we define a patient and the actions he will take, like requesting beds, start or finishing the internment or giving up if there is no bed:

def patient(env, name, lt):
    global used_beds
    global days    
    with lt.bed.request() as request:
        results = yield request | env.timeout(0)

        if request in results:
            print("%s interns at day %s" % (name, int(env.now)))
            yield env.process(lt.internment(name))
            print("%s leaves the bed at day %s" % (name, int(env.now)))
        else:
            print("%s had no bed at day %s" % (name, int(env.now)))

    used_beds.append(lt.bed.count)
    days.append(int(env.now))

The 'used_beds' and 'days' variables saves the bed occupancy and respective days. The following step is to define the numbers of patients that request system care for each day. These values we get from the infected curve I(t) considering 20% of the infecteds, subtracting the number of previous day, defining new cases as 20% of I(t) - I(t-1). When the curve starts do decline there are no new cases and we consider new cases = 0, otherwise the number will be negative.

def setup(env, num_beds, intern_time):
    bed = Healthcare(env, num_beds, intern_time)
    j = 0
    for t in range(130):
        global intern_total
        new_cases = int((I[t] - I[t - 1]) * 0.2)
        if new_cases < 0:
            new_cases = 0
        intern_total += new_cases
        for i in range(new_cases):
            j += 1
            env.process(patient(env, 'Patient %d' % j, bed))
        yield env.timeout(1)

Now, it is time to run the simulation:

random.seed(RANDOM_SEED)
env = simpy.Environment()
env.process(setup(env, NUM_BEDS, INTERN_TIME))
env.run(until=SIM_TIME)

And finally, plot the graph with the curves, choosing the curve we want by comment or uncomment the respective line of plotting.

# Plot data
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)

# The 3 curves S(t), I(t) and R(t)
# ax.plot(t, S, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I, 'm', alpha=0.5, lw=2, label='Infected')
ax.plot(t, I * 0.2, 'r', alpha=0.5, lw=2, label='20% of infected need hospital beds')
# ax.plot(t, R, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
plt.axhline(y=NUM_BEDS, color='b', alpha=0.5, lw=2, label='number of available beds')
ax.plot(occup, alpha=0.5, color='g', lw=2, label='bed occupancy')

ax.set_xlabel('days')
ax.set_ylabel('people')
plt.title('Epidemic curve')
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
plt.ticklabel_format(axis="y", style="plain")
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()
sys.stdout.close()

See the resulting curve:

Notice the green curve. It is the bed occupancy and it does not reach the number of available bed's line. It means that, although the red curve crosses the bed's line, in fact it is possible to attend the needed patients considering that some of the beds will be available as the patients are being recovered from treatment, so this give us a more realistic bed occupancy and requirement and help us to better planning the resources needed. Try to change the period of internment with values different from 15 and 21 and see the results... It is a good exercise.

Take a look also in the txt file saved in the folder. You will see the dynamic of the simulation, something like that:

...
Patient 251 interns at day 31
Patient 252 interns at day 31
Patient 253 interns at day 31
Patient 254 interns at day 31
Patient 5 leaves the bed at day 31
Patient 6 leaves the bed at day 31
Patient 255 interns at day 32
Patient 256 interns at day 32

This calculation could help us to better plan the resources needed to face the pandemia and it is also possible to add more resources like respirators and also human resources taking more information into consideration. It is a helpful tool!

Joining all the code, we have:

# ******************************************************************************
# * Sizing the healthcare system resources to face the pandemia
# * Author: Sergio Bastos    04.17.2020
# ******************************************************************************

import numpy as np
import scipy
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import pandas as pd
import sys
import random
import simpy

# historical data
hist = pd.Series([5, 5, 5, 5, 7, 9, 14, 17, 27, 29, 31, 40, 46, 56, 66, 71])

# Total population N, gamma
N = 1000000
gamma = 0.5

# Exponential adjusting of the historical data
popt_exponential, pcov_exponential = \
    scipy.optimize.curve_fit(lambda t, a, b: a * np.exp(b * t),
                             hist.index, hist, p0=(5, 0.3))

# Parameters for the SIR model
I0 = popt_exponential[0]
beta = gamma + popt_exponential[1]
R0 = beta / gamma
S0 = N - I0 - R0

t = np.linspace(0, 130, 130)


# Diferential equations of the SIR model
def deriv(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt


# Initial conditions vector
y0 = S0, I0, R0

# Integration
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T

# Simulation

RANDOM_SEED = 100
NUM_BEDS = 7000
INTERN_TIME = random.triangular(15, 21)
SIM_TIME = 140

sys.stdout = open("./sim_dynamics.txt", "w")


class Healthcare(object):
    def __init__(self, env, num_beds, intern_time):
        self.env = env
        self.bed = simpy.Resource(env, num_beds)
        self.intern_time = intern_time

    def internment(self, patient):
        yield self.env.timeout(INTERN_TIME)


def patient(env, name, lt):
    global used_beds
    global days
    with lt.bed.request() as request:
        results = yield request | env.timeout(0)

        if request in results:
            print("%s interns at day %s" % (name, int(env.now)))
            yield env.process(lt.internment(name))
            print("%s leaves the bed at day %s" % (name, int(env.now)))
        else:
            print("%s stays without care at day %s" % (name, int(env.now)))

    used_beds.append(lt.bed.count)
    days.append(int(env.now))


def setup(env, num_beds, intern_time):
    bed = Healthcare(env, num_beds, intern_time)
    j = 0
    for t in range(130):
        global intern_total
        new_cases = int((I[t] - I[t - 1]) * 0.2)
        if new_cases < 0:
            new_cases = 0
        intern_total += new_cases
        for i in range(new_cases):
            j += 1
            env.process(patient(env, 'Patient %d' % j, bed))
        yield env.timeout(1)


intern_total = 0
used_beds = []
days = []

random.seed(RANDOM_SEED)
env = simpy.Environment()
env.process(setup(env, NUM_BEDS, INTERN_TIME))
env.run(until=SIM_TIME)

occup = pd.Series(used_beds, index=days)
print(occup)
occup = occup.loc[~occup.index.duplicated(keep='first')]
print(occup)

# Plot data
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)

# The 3 curves S(t), I(t) and R(t)
# ax.plot(t, S, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I, 'm', alpha=0.5, lw=2, label='Infected')
ax.plot(t, I * 0.2, 'r', alpha=0.5, lw=2, label='20% of infected need hospital beds')
# ax.plot(t, R, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
plt.axhline(y=NUM_BEDS, color='b', alpha=0.5, lw=2, label='number of available beds')
ax.plot(occup, alpha=0.5, color='g', lw=2, label='bed occupancy')

ax.set_xlabel('days')
ax.set_ylabel('people')
plt.title('Epidemic curve')
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
plt.ticklabel_format(axis="y", style="plain")
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()
sys.stdout.close()

References:

[1] W. O. Kermack e A. G. McKendrick, editor. A Contribution to The Mathematical Theory of Epidemics, volume 115 of A. Royal Society of London, August 1927.

[2] Liangrong Peng, Wuyue Yang, Dongyan Zhang, Changjing Zhuge, and Liu Hong. Epidemic analysis of covid-19 in china by dynamical modeling. https://www.medrxiv.org/content/10.1101/2020.02.16.20023465v1.full.pdf, 2020.

[3] Christian Hill. Learning Scientific Programming with Python. Cambridge UniversityPress, 2016.

[4] SimPy documentation https://simpy.readthedocs.io/en/latest/contents.html