|
| 1 | +import pandas as pd |
| 2 | +import streamlit as st |
| 3 | +import numpy as np |
| 4 | +import matplotlib |
| 5 | +matplotlib.use('Agg') |
| 6 | +import matplotlib.pyplot as plt |
| 7 | + |
| 8 | +hide_menu_style = """ |
| 9 | + <style> |
| 10 | + #MainMenu {visibility: hidden;} |
| 11 | + </style> |
| 12 | + """ |
| 13 | +st.markdown(hide_menu_style, unsafe_allow_html=True) |
| 14 | + |
| 15 | +delaware = 564696 |
| 16 | +chester = 519293 |
| 17 | +montgomery = 826075 |
| 18 | +bucks = 628341 |
| 19 | +philly = 1581000 |
| 20 | +S_default = delaware + chester + montgomery + bucks + philly |
| 21 | +known_infections = 31 |
| 22 | + |
| 23 | +# Widgets |
| 24 | +initial_infections = st.sidebar.number_input( |
| 25 | + "Currently Known Regional Infections", value=known_infections, step=10, format="%i" |
| 26 | +) |
| 27 | +current_hosp = st.sidebar.number_input( |
| 28 | + "Currently Hospitalized COVID-19 Patients", value=2, step=1, format="%i" |
| 29 | +) |
| 30 | +doubling_time = st.sidebar.number_input( |
| 31 | + "Doubling Time (days)", value=6, step=1, format="%i" |
| 32 | +) |
| 33 | +hosp_rate = ( |
| 34 | + st.sidebar.number_input("Hospitalization %", 0, 100, value=5, step=1, format="%i") |
| 35 | + / 100.0 |
| 36 | +) |
| 37 | +icu_rate = ( |
| 38 | + st.sidebar.number_input("ICU %", 0, 100, value=2, step=1, format="%i") / 100.0 |
| 39 | +) |
| 40 | +vent_rate = ( |
| 41 | + st.sidebar.number_input("Ventilated %", 0, 100, value=1, step=1, format="%i") |
| 42 | + / 100.0 |
| 43 | +) |
| 44 | +hosp_los = st.sidebar.number_input("Hospital LOS", value=7, step=1, format="%i") |
| 45 | +icu_los = st.sidebar.number_input("ICU LOS", value=9, step=1, format="%i") |
| 46 | +vent_los = st.sidebar.number_input("Vent LOS", value=10, step=1, format="%i") |
| 47 | +Penn_market_share = ( |
| 48 | + st.sidebar.number_input( |
| 49 | + "Hospital Market Share (%)", 0.0, 100.0, value=15.0, step=1.0, format="%f" |
| 50 | + ) |
| 51 | + / 100.0 |
| 52 | +) |
| 53 | +S = st.sidebar.number_input( |
| 54 | + "Regional Population", value=S_default, step=100000, format="%i" |
| 55 | +) |
| 56 | + |
| 57 | +total_infections = current_hosp / Penn_market_share / hosp_rate |
| 58 | +detection_prob = initial_infections / total_infections |
| 59 | + |
| 60 | +st.title("COVID-19 Hospital Impact Model for Epidemics") |
| 61 | +st.markdown( |
| 62 | + """*This tool was developed by the [Predictive Healthcare team](http://predictivehealthcare.pennmedicine.org/) at Penn Medicine. For questions and comments please see our [contact page](http://predictivehealthcare.pennmedicine.org/contact/).*""" |
| 63 | +) |
| 64 | + |
| 65 | +if st.checkbox("Show more info about this tool"): |
| 66 | + st.subheader( |
| 67 | + "[Discrete-time SIR modeling](https://mathworld.wolfram.com/SIRModel.html) of infections/recovery" |
| 68 | + ) |
| 69 | + st.markdown( |
| 70 | + """The model consists of individuals who are either _Susceptible_ ($S$), _Infected_ ($I$), or _Recovered_ ($R$). |
| 71 | +
|
| 72 | +The epidemic proceeds via a growth and decline process. This is the core model of infectious disease spread and has been in use in epidemiology for many years.""" |
| 73 | + ) |
| 74 | + st.markdown("""The dynamics are given by the following 3 equations.""") |
| 75 | + |
| 76 | + st.latex("S_{t+1} = (-\\beta S_t I_t) + S_t") |
| 77 | + st.latex("I_{t+1} = (\\beta S_t I_t - \\gamma I_t) + I_t") |
| 78 | + st.latex("R_{t+1} = (\\gamma I_t) + R_t") |
| 79 | + |
| 80 | + st.markdown( |
| 81 | + """To project the expected impact to Penn Medicine, we estimate the terms of the model. |
| 82 | +
|
| 83 | +To do this, we use a combination of estimates from other locations, informed estimates based on logical reasoning, and best guesses from the American Hospital Association. |
| 84 | +
|
| 85 | +
|
| 86 | +### Parameters |
| 87 | +First, we need to express the two parameters $\\beta$ and $\\gamma$ in terms of quantities we can estimate. |
| 88 | +
|
| 89 | +- The $\\gamma$ parameter represents 1 over the mean recovery time in days. Since the CDC is recommending 14 days of self-quarantine, we'll use $\\gamma = 1/14$. |
| 90 | +- Next, the AHA says to expect a doubling time $T_d$ of 7-10 days. That means an early-phase rate of growth can be computed by using the doubling time formula: |
| 91 | +""" |
| 92 | + ) |
| 93 | + st.latex("g = 2^{1/T_d} - 1") |
| 94 | + |
| 95 | + st.markdown( |
| 96 | + """ |
| 97 | +- Since the rate of new infections in the SIR model is $g = \\beta S - \\gamma$, and we've already computed $\\gamma$, $\\beta$ becomes a function of the initial population size of susceptible individuals. |
| 98 | +$$\\beta = (g + \\gamma)/s$$ |
| 99 | +
|
| 100 | +### Initial Conditions |
| 101 | +
|
| 102 | +- The total size of the susceptible population will be the entire catchment area for Penn Medicine entities (HUP, PAH, PMC, CCH) |
| 103 | + - Delaware = {delaware} |
| 104 | + - Chester = {chester} |
| 105 | + - Montgomery = {montgomery} |
| 106 | + - Bucks = {bucks} |
| 107 | + - Philly = {philly} |
| 108 | +- The initial number of infected will be the total number of confirmed cases in the area ({initial_infections}), divided by some detection probability to account for under testing {detection_prob:.2f}.""".format( |
| 109 | + delaware=delaware, |
| 110 | + chester=chester, |
| 111 | + montgomery=montgomery, |
| 112 | + bucks=bucks, |
| 113 | + philly=philly, |
| 114 | + initial_infections=initial_infections, |
| 115 | + detection_prob=detection_prob, |
| 116 | + ) |
| 117 | + ) |
| 118 | + |
| 119 | +# The SIR model, one time step |
| 120 | +def sir(y, beta, gamma, N): |
| 121 | + S, I, R = y |
| 122 | + Sn = (-beta * S * I) + S |
| 123 | + In = (beta * S * I - gamma * I) + I |
| 124 | + Rn = gamma * I + R |
| 125 | + if Sn < 0: |
| 126 | + Sn = 0 |
| 127 | + if In < 0: |
| 128 | + In = 0 |
| 129 | + if Rn < 0: |
| 130 | + Rn = 0 |
| 131 | + |
| 132 | + scale = N / (Sn + In + Rn) |
| 133 | + return Sn * scale, In * scale, Rn * scale |
| 134 | + |
| 135 | + |
| 136 | +# Run the SIR model forward in time |
| 137 | +def sim_sir(S, I, R, beta, gamma, n_days, beta_decay=None): |
| 138 | + N = S + I + R |
| 139 | + s, i, r = [S], [I], [R] |
| 140 | + for day in range(n_days): |
| 141 | + y = S, I, R |
| 142 | + S, I, R = sir(y, beta, gamma, N) |
| 143 | + if beta_decay: |
| 144 | + beta = beta * (1 - beta_decay) |
| 145 | + s.append(S) |
| 146 | + i.append(I) |
| 147 | + r.append(R) |
| 148 | + |
| 149 | + s, i, r = np.array(s), np.array(i), np.array(r) |
| 150 | + return s, i, r |
| 151 | + |
| 152 | + |
| 153 | +## RUN THE MODEL |
| 154 | + |
| 155 | +S, I, R = S, initial_infections / detection_prob, 0 |
| 156 | + |
| 157 | +intrinsic_growth_rate = 2 ** (1 / doubling_time) - 1 |
| 158 | + |
| 159 | +recovery_days = 14.0 |
| 160 | +# mean recovery rate, gamma, (in 1/days). |
| 161 | +gamma = 1 / recovery_days |
| 162 | + |
| 163 | +# Contact rate, beta |
| 164 | +beta = ( |
| 165 | + intrinsic_growth_rate + gamma |
| 166 | +) / S # {rate based on doubling time} / {initial S} |
| 167 | + |
| 168 | + |
| 169 | +n_days = st.slider("Number of days to project", 30, 200, 60, 1, "%i") |
| 170 | + |
| 171 | +beta_decay = 0.0 |
| 172 | +s, i, r = sim_sir(S, I, R, beta, gamma, n_days, beta_decay=beta_decay) |
| 173 | + |
| 174 | + |
| 175 | +hosp = i * hosp_rate * Penn_market_share |
| 176 | +icu = i * icu_rate * Penn_market_share |
| 177 | +vent = i * vent_rate * Penn_market_share |
| 178 | + |
| 179 | +days = np.array(range(0, n_days + 1)) |
| 180 | +data_list = [days, hosp, icu, vent] |
| 181 | +data_dict = dict(zip(["day", "hosp", "icu", "vent"], data_list)) |
| 182 | + |
| 183 | +projection = pd.DataFrame.from_dict(data_dict) |
| 184 | + |
| 185 | +st.subheader("New Admissions") |
| 186 | +st.markdown("Projected number of **daily** COVID-19 admissions at Penn hospitals") |
| 187 | + |
| 188 | +# New cases |
| 189 | +projection_admits = projection.iloc[:-1, :] - projection.shift(1) |
| 190 | +projection_admits[projection_admits < 0] = 0 |
| 191 | + |
| 192 | +plot_projection_days = n_days - 10 |
| 193 | +projection_admits["day"] = range(projection_admits.shape[0]) |
| 194 | + |
| 195 | +fig, ax = plt.subplots(1, 1, figsize=(10, 4)) |
| 196 | +ax.plot( |
| 197 | + projection_admits.head(plot_projection_days)["hosp"], ".-", label="Hospitalized" |
| 198 | +) |
| 199 | +ax.plot(projection_admits.head(plot_projection_days)["icu"], ".-", label="ICU") |
| 200 | +ax.plot(projection_admits.head(plot_projection_days)["vent"], ".-", label="Ventilated") |
| 201 | +ax.legend(loc=0) |
| 202 | +ax.set_xlabel("Days from today") |
| 203 | +ax.grid("on") |
| 204 | +ax.set_ylabel("Daily Admissions") |
| 205 | +st.pyplot() |
| 206 | + |
| 207 | +admits_table = projection_admits[np.mod(projection_admits.index, 7) == 0].copy() |
| 208 | +admits_table["day"] = admits_table.index |
| 209 | +admits_table.index = range(admits_table.shape[0]) |
| 210 | +admits_table = admits_table.fillna(0).astype(int) |
| 211 | + |
| 212 | +if st.checkbox("Show Projected Admissions in tabular form"): |
| 213 | + st.dataframe(admits_table) |
| 214 | + |
| 215 | +st.subheader("Admitted Patients (Census)") |
| 216 | +st.markdown( |
| 217 | + "Projected **census** of COVID-19 patients, accounting for arrivals and discharges at Penn hospitals" |
| 218 | +) |
| 219 | + |
| 220 | +# ALOS for each category of COVID-19 case (total guesses) |
| 221 | + |
| 222 | +los_dict = { |
| 223 | + "hosp": hosp_los, |
| 224 | + "icu": icu_los, |
| 225 | + "vent": vent_los, |
| 226 | +} |
| 227 | + |
| 228 | +fig, ax = plt.subplots(1, 1, figsize=(10, 4)) |
| 229 | + |
| 230 | +census_dict = {} |
| 231 | +for k, los in los_dict.items(): |
| 232 | + census = ( |
| 233 | + projection_admits.cumsum().iloc[:-los, :] |
| 234 | + - projection_admits.cumsum().shift(los).fillna(0) |
| 235 | + ).apply(np.ceil) |
| 236 | + census_dict[k] = census[k] |
| 237 | + ax.plot(census.head(plot_projection_days)[k], ".-", label=k + " census") |
| 238 | + ax.legend(loc=0) |
| 239 | + |
| 240 | +ax.set_xlabel("Days from today") |
| 241 | +ax.grid("on") |
| 242 | +ax.set_ylabel("Census") |
| 243 | +st.pyplot() |
| 244 | + |
| 245 | +census_df = pd.DataFrame(census_dict) |
| 246 | +census_df["day"] = census_df.index |
| 247 | +census_df = census_df[["day", "hosp", "icu", "vent"]] |
| 248 | + |
| 249 | +census_table = census_df[np.mod(census_df.index, 7) == 0].copy() |
| 250 | +census_table.index = range(census_table.shape[0]) |
| 251 | +census_table.loc[0, :] = 0 |
| 252 | +census_table = census_table.dropna().astype(int) |
| 253 | + |
| 254 | +if st.checkbox("Show Projected Census in tabular form"): |
| 255 | + st.dataframe(census_table) |
| 256 | + |
| 257 | +st.markdown( |
| 258 | + """**Click the checkbox below to view additional data generated by this simulation**""" |
| 259 | +) |
| 260 | +if st.checkbox("Show Additional Projections"): |
| 261 | + st.subheader( |
| 262 | + "The number of infected and recovered individuals in the hospital catchment region at any given moment" |
| 263 | + ) |
| 264 | + fig, ax = plt.subplots(1, 1, figsize=(10, 4)) |
| 265 | + ax.plot(i, label="Infected") |
| 266 | + ax.plot(r, label="Recovered") |
| 267 | + ax.legend(loc=0) |
| 268 | + ax.set_xlabel("days from today") |
| 269 | + ax.set_ylabel("Case Volume") |
| 270 | + ax.grid("on") |
| 271 | + st.pyplot() |
| 272 | + |
| 273 | + # Show data |
| 274 | + days = np.array(range(0, n_days + 1)) |
| 275 | + data_list = [days, s, i, r] |
| 276 | + data_dict = dict(zip(["day", "susceptible", "infections", "recovered"], data_list)) |
| 277 | + projection_area = pd.DataFrame.from_dict(data_dict) |
| 278 | + infect_table = (projection_area.iloc[::7, :]).apply(np.floor) |
| 279 | + infect_table.index = range(infect_table.shape[0]) |
| 280 | + |
| 281 | + if st.checkbox("Show Raw SIR Similation Data"): |
| 282 | + st.dataframe(infect_table) |
| 283 | + |
| 284 | +st.subheader("References & Acknowledgements") |
| 285 | +st.markdown( |
| 286 | + """* AHA Webinar, Feb 26, James Lawler, MD, an associate professor University of Nebraska Medical Center, What Healthcare Leaders Need To Know: Preparing for the COVID-19 |
| 287 | +* We would like to recognize the valuable assistance in consultation and review of model assumptions by Michael Z. Levy, PhD, Associate Professor of Epidemiology, Department of Biostatistics, Epidemiology and Informatics at the Perelman School of Medicine |
| 288 | + """ |
| 289 | +) |
| 290 | +st.markdown("© 2020, The Trustees of the University of Pennsylvania") |
0 commit comments