10 Hospital systems play a critical role in treating injuries and preventing additional deaths 11 during disaster emergency response. Natural disasters hinder the ability of hospital systems to 12 operate at full capacity. Therefore, it is important for cities to develop policies and standards 13 that enable hospitals' continuous operations to provide patients with timely treatment and en- 14 sure urban resilience. Here, we present a methodology to evaluate emergency response based 15 on a probabilistic model that assesses the loss of hospital functions and quanti es multiseverity 16 injuries as a result of earthquake damage. The proposed methodology is able to design e ective 17 plans for patient transferal and allocation of medical resources using an optimization formu- 18 lation. This methodology is applied to Lima, Peru, subjected to a disaster scenario based on 19 the M 8.0 earthquake that occurred there in 1940. Our results show that the spatial distribu- 20 tion of health service demands mismatches the post-earthquake capacities of hospitals, leaving 21 large zones on the periphery of Lima signi cantly underserved. This study demonstrates how 22 emergency plans that leverage hospital-system coordination can address this demand-capacity 23 mismatch, enabling e ective patient transfers, ambulance usage, and deployment of emergency 24 medical teams. 25 Hospital systems are at the core of disaster resilience because they must provide timely critical 26 healthcare services to communities during and after an emergency response.1Because cities are 27 becoming larger and more densely populated, natural disasters are impacting public health on a 28 larger scale. A database including the most 21,000 devastating disasters worldwide since 1900 29 indicates that 50% of disasters with the largest number of injuries occurred only during the last 20 30 years.2Natural disasters such as earthquakes, landslides, oods, typhoons put heavy demands on 31 hospital systems because these disasters can cause thousands or even tens of thousands of injuries 32 in a short timespan (Figure 1). At the same time, natural disasters cause massive disruptions to 33 hospital systems by damaging their supporting infrastructure. For example, the M 7.6 1999 Turkey 34 Corresponding author: Luis Ceferino. E-mail: lceferinor@gmail.com Preprint submitted to Nature Communications 1 earthquake caused around 50,000 injuries in Izmit and disrupted 10 major hospitals, which required 35 relocation of most patients from these hospitals.336 Because hospital systems are so critical, the World Health Organization (WHO) and Pan- 37 American Health Organization (PAHO) urge countries to institute policies to strengthen hospi- 38 tal capacities and enhance coordination in the hospital system to make ecient use of resources 39 at national and regional levels during emergency response.4,5To e ectively develop measures 40 for capacity-enhancing prioritization and resource sharing and allocation, national and regional 41 governments require information based on robust methodologies that can characterize hospitals' 42 emergency response as an interconnected system on a large urban scale. 43 However, most previous studies have primarily focused on modeling emergency response only at 44 single-hospital scale as opposed to characterizing the response of hospital systems on a large urban 45 scale. Some of these studies relied on disaster analytics to evaluate post-disaster functionality of the 46 supporting infrastructure in the individual hospitals.6{9Other studies used emergency medicine 47 modeling tools, such as discrete event simulation (DES) and ow models, to characterize emergency 48 response and evaluate post-disaster resource allocation but also at a single-hospital scale.10{13Lack 49 of methods and high-resolution disaster risk data have hindered the extension of single-hospital 50 scale analyses to system-level analyses on an urban scale. As a result, regional emergency response 51 policies have not e ectively addressed capacity-enhancing prioritization and resource sharing and 52 allocation in hospital systems, especially in large and complex urban centers. 53 Here, we present ndings from a methodology that characterizes the disaster emergency response 54 of hospital systems on a large urban scale. Our integrative methodology combines models of 55 multiseverity earthquake casualty estimation14,15and post-earthquake hospital functionality with a 56 proposed network ow model for hospital systems. We focus on seismic hazard because earthquakes 57 are the natural disasters that have caused the largest number of injuries in most countries (Figure 1). 58 The methodology is applied to Lima, Peru, based on a M 8.0 earthquake and includes an evaluation 59 of e ective emergency plans for allocation of hospital resources and patient transfers. We selected 60 Lima because it has a high seismic risk and it has recently built a unique dataset containing high- 61 resolution hospital vulnerability. We use citywide data on the seismic vulnerability of more than 62 1.5 M buildings in Lima to estimate casualties and data including the seismic vulnerability of 41 63 public hospital campuses (composed of +700 buildings) and their respective operating rooms and 64 ambulance resources.9,16,1765 We propose a metric based on patient waiting times and e ective use of ambulance patient 66 transfers as a performance measure for developing emergency response plans. Our focus is on 67 high-severity injuries that require surgical procedures. We evaluate the spatial distributions of 68 high-severity injuries in the city at a higher spatial resolution (i.e., 1kmx1km) than other widely 69 used methods.18Then, we compare the spatial distribution of casualties with the distribution of 70 functional operating rooms in the hospital system, identifying the zones more likely to be under- 71 served during the emergency response. Combining the network ow model with an optimization 72 formulation, we assess the performance of four alternative emergency response plans to treat the 73 patients in the city. 74 The rst and second emergency plans are baseline strategies with limited levels of coordination. 75 In both strategies, hospitals without available operating resources will use their own ambulances to 76 transfer patients, but in the rst strategy, patients will only be sent to the closest working hospital, 77 whereas in the second strategy, patients will be sent to the hospital with the largest number of 78 functional operating rooms. 79 2 The third and fourth emergency plans are strategies with higher levels of coordination that use 80 the optimization formulation on the system performance metric. In both strategies, all hospitals 81 will share ambulance resources across the system to transfer patients according to post-earthquake 82 needs, but in the third strategy, the system only uses the residual operating room capacities in 83 the hospitals, whereas in the fourth strategy, emergency medical teams (EMTs) supply the system 84 with additional mobile operating rooms in key locations in the city. Through the Action Plan for 85 Humanitarian Assistance, the WHO and PAHO require countries to elaborate policies for deploying 86 ETMs to assist people a ected by emergencies and disasters,19,20thus, this study aims to directly 87 inform policies for EMT deployment in countries with high seismic risk. 88 We analyze the behavior and performance of these four plans during the emergency response and 89 discuss their implications in terms of patient treatment times, ambulance usage, and patient trans- 90 fers. We utilize trac data to show the most important roads for patient transfers from localized 91 zones with lower hospital capacity to zones with higher capacity in the city. This research repre- 92 sents a rst-cut assessment on the e ectiveness of emergency response policies to inform city-scale 93 decision-making that leads to more e ective treatment of patients during an emergency response 94 to a major earthquake. 95 Figure 1: Per-country distribution of disasters with the highest number of injuries since 1900. The sizes of the circles indicate the relative number of injuries and the colors indicate the natural disaster type. Since 1900, 117 countries that experienced at least one natural disaster with more than 100 injuries. Earthquakes were the natural disaster with largest injury tolls in 45 of these countries, followed by storms in 38 of these countries. Earthquake events can cause large number of injuries suddenly. For example, the 2008 Great Sichuan earthquake injured more than 368k people in only two minutes. Data from EM-DAT.2 1 Results 96 We applied our methodology to Lima, a large city with a population close to 10 million people,2197 where previous large earthquakes have caused large numbers of casualties.22,23Because the last 98 3 large earthquake in the region occurred more than 40 years ago, studies indicate that the city is 99 currently exposed to high hazard of large-magnitude earthquakes.24,25Using our proposed method- 100 ology, we characterize the emergency response and evaluate response plans to e ectively to treat 101 patients in Lima after an earthquake scenario of large magnitude occurring in the near future. This 102 earthquake scenario was simulated according to the seismotectonics of the 1940 M 8.0 earthquake, 103 which occurred in close proximity to Lima.26Figure 2 shows the estimated rupture area of the 1940 104 earthquake and its proximity to the city. Our methodology estimates the impact of this disaster 105 scenario on the demands on healthcare by quantifying earthquake casualties and on the capacity 106 of healthcare by quantifying the post-earthquake reduction in functionality in the hospital system. 107 The results of applying our methodology are discussed here and the details of the methodology 108 formulation, work ow and required data are described in the Methods section. 109 Figure 2: Earthquake scenario representing the M 8.0 1940 earthquake in Lima. The earthquake occurred in the subduction fault in the coast of Lima and caused widespread damage to the city.22, 27The estimated area of fault rupture is shown in red. The edge dimensions were estimated with empirical formulas.28The fault plane dips 15, where the edge underneath the coast is deeper than the edge under the ocean. The median peak ground acceleration (PGA) is also estimated with empirical formulas.29The shaking attenuates for regions further away from the rupture in the fault plane. Lima city and its districts are delimited by the black shapes. 1.1 Earthquake Casualties 110 We found that on average close to 4.7k people will require surgical procedures in operating rooms 111 after the M 8.0 earthquake. This estimate results from applying a probabilistic model that utilizes 112 high-resolution building seismic vulnerability data, population distribution and soil conditions to 113 4 evaluate multiseverity earthquake casualties caused by widespread building damage14,15(see Meth- 114 ods). Earthquake injuries can have di erent severity degrees, ranging from small bruises to more 115 serious spinal cord injuries.30{34The 4.7k patients requiring surgical procedures will have high 116 severity injuries such us compound bone fractures, punctured organs or crush syndrome with open 117 wounds, thus they require timely interventions for stabilization and treatment. 118 Our results are designed for a nighttime scenario, when most people are inside residential build- 119 ings, because residential infrastructure is particularly vulnerable in Lima. Predominantly the city's 120 periphery has vulnerable residential infrastructure as a result of poor construction practices and 121 lack of seismic code enforcement.35,36Figure 3a shows the spatial distribution of the average num- 122 ber of patients that will require the surgical procedures. A comparison with the spatial distribution 123 of nighttime population density in Lima (Figure 9) indicates that most of these patients are located 124 in high-density zones. However, the ratio between the number of injured people and the total num- 125 ber of people follows a di erent pattern (Figure 3b). The spatial distribution of this ratio re ects 126 the uneven distribution of ground shaking intensities and seismic vulnerabilities of buildings in the 127 city. It shows that people living closer to the coastline and in the city's peripheral zones have 128 higher earthquake injury risk as a result of higher ground shaking (Figure 2) and more vulnerable 129 buildings, respectively. 130 (a) (b) Figure 3: Casualty scenario for M 8.0 earthquake occurring at nighttime in Lima. (a) Spatial distribution in km2of earthquake injuries requiring surgical procedures after the M 8.0 seismic event. (b) Spatial distribution of earthquake injury ratios, i.e., number of injuries as a percentage of the population per km2. The intervals in the two plots represent quintiles (5-quantile) on the spatial data. 5 1.2 Post-earthquake Hospital Capacity 131 We found that on average only 93 of 182 total hospital operating rooms (51%) will be functional after 132 the M 8.0 earthquake. This estimate results from performing a probabilistic earthquake simulation 133 on a high-resolution dataset (see Methods). The dataset includes the structural vulnerabilities of 134 +700 buildings belonging to 41 healthcare campuses,16the operating room resources, and \Hospital 135 Safety Index" (HSI) of each campus. HSI is a metric created by WHO to measure post-disaster 136 functionality potential due to multiple factors such as backup water, power, medical resources and 137 hospital accessibility.37This unique dataset in combination with the earthquake simulation enables 138 us to capture the residual hospital functionality on a large urban scale. 139 Figure 5a shows the spatial distribution of both the operating rooms in the dataset and the 140 average predictions of operating rooms after the earthquake. Both spatial distributions are heavily 141 uneven across the city. In the dataset, 95 operating rooms (52%) are concentrated in only four 142 centric districts, Lima, Bre~ na, La Victoria and Jes us Mar a, whose summed areas represent less 143 than 2% than the total area of the city. The earthquake slightly worsens such a resource central- 144 ization due to the non-uniform spatial distribution of earthquake shaking and the variations in the 145 vulnerabilities of hospitals' buildings according to their construction age and standards or struc- 146 tural types (Figure 11). As a result, we estimate that these four districts will have 55 functional 147 operating rooms, 59% of the total functional operating rooms, in the emergency response. 148 Additionally, addressing the centralization issue can become even more critical because the 149 number of injuries needing surgical procedures and the functional operating rooms are negatively 150 correlated. Our ndings show a strong correlation (-0.49) between the simulations of earthquake 151 injuries and functional operating rooms across the system in the city (Figure 4). Such a large 152 correlation indicates that an earthquake that injures a larger amount of people will likely be very 153 destructive; thus, that scenario will also cause a heavier disruption to the hospital system. 154 1.3 Demand-capacity Mismatch of Health Services 155 We analyzed patient arrivals to the hospitals and found that patient distribution signi cantly 156 mismatches the distribution of residual hospital resources after the earthquake. We assumed that 157 search and rescue (SAR) teams, relatives, friends and neighbors will initially transport patients 158 to the triage areas in the closest hospitals as it occurred after previous earthquakes.38Figure 5b 159 shows the distribution of the total arrivals of the patients who will need surgical procedures in each 160 hospital. In contrast to the distribution of functional operating rooms, injuries are mainly located 161 in the periphery. Only 596 patients would arrive at the hospitals in the four centric districts 162 highlighted in Figure 5b. Those patients represent only 13% of the total demand for surgical 163 procedures. However, as described earlier, these four districts will concentrate 59% of the functional 164 operating rooms available. 165 Such a mismatch in the distribution of demands and capacities creates localized health service 166 imbalances leading to long patient waiting times, with particularly severe e ects in the periphery. 167 Emergency plans can play a key role in addressing this mismatch and improve the treatment 168 e ectiveness in the city if they either mobilize patients from lower-capacity zones to higher-capacity 169 zones or supply lower-capacity zones with additional resources. We tested four emergency plans. 170 Two of them are baseline strategies that only require limited coordination in the system, whereas 171 the other two are strategies that require higher coordination at the system level. To evaluate their 172 performance during the emergency response, we used a system metric based on city-wide patient 173 6 Figure 4: 1,000 simulations of number of casualties needing surgical procedures in operating rooms (ORs) and number of functional operating rooms after the earthquake for the M 8.0 earthquake. The simulations result from probabilistic earthquake modeling (see Methods) and capture uncertainty in ground shaking, building damage, injury occurrence and hospital functionality. The linear trend indicates a negative corre- lation between the functional ORs and the number of injuries in the simulations. waiting times and e ective use of ambulance resources (see Methods). Waiting times are a key 174 metric to establish necessary patient stabilization procedures until there is an available operating 175 room in the queue, and the use of ambulances is a complementary measure to ensure that patient 176 transfers occur e ectively. Combining simulations on post-earthquake demand-capacity with data 177 including the available ambulances at each campus and regular trac conditions in Lima from 178 Google Maps API, we conducted a probabilistic evaluation of the system metric performance for 179 the four emergency plans. 180 1.4 Baseline Strategies with Limited Coordination 181 In the rst strategy, hospitals send patients to the closest hospital with functional operating rooms 182 only if all their operating rooms are non-functional after the earthquake. In this strategy, hospitals 183 use their own ambulance resources to transfer their patients. This strategy only requires limited 184 coordination between pairs of hospitals located relatively close to each other, representing an emer- 185 gency response where the system becomes a set of islands composed of districts or neighborhoods 186 that treat injuries independently of each other. With this strategy, our mean estimates indicate 187 that the average waiting time will be 30 days to receive treatment in operating rooms (Figure 188 6). In some worst-case scenarios, this metric could even increase up to 76 days, as indicated by 189 the 90th-percentile of this distribution. The 1988 Armenia39and the 1991 Costa Rica38earth- 190 7 (a) (b) Figure 5: Distribution of operating rooms in Lima. The circle size represents the relative number of operating rooms. (a) Current number of operating rooms in hospital locations16and mean estimates of functional operating rooms after the M 8.0 earthquake. (b) Mean estimates of total arrivals of patients who will need surgical procedures after the earthquake. quakes showed that delayed surgical treatment can worsen the patients' health to life-threatening 191 conditions, for example, those who need resuscitative surgery, e.g., intra-abdominal hemorrhage 192 or emergency amputation. Thus, such long waiting times in Lima can result in many additional 193 deaths. 194 For the rst strategy, Figure 7a shows the mean estimates of the spatial distribution of treated 195 patients at each hospital and the patient transfers between hospital pairs. Though this strategy can 196 ooad demands in critical zones, it does not e ectively mobilize patients from the lower-capacity 197 zones to higher-capacity zones. Hospitals with more functional operating rooms treat a similar 198 number of patients as hospitals with fewer operating rooms. Hospitals in the four centric districts 199 highlighted previously only treat 1.3k patients, which represents 28% of the demand for operating 200 rooms, despite having 59% of the total capacity. Additionally, because hospitals do not share 201 ambulance resources, we nd that ambulances are the bottleneck of the system in the periphery. 202 Hospitals with limited ambulances have to transfer large numbers of patients to ooad the high 203 demand for operating rooms. Thus, if they do not have sucient ambulance resources, their patients 204 will lose the opportunity to be treated more promptly in other less crowded hospitals. 205 In the second strategy, hospitals send patients to the hospital with the largest number of func- 206 tional operating rooms. In this strategy, hospitals send patients only if all their operating rooms are 207 non-functional and use their own ambulance resources. With this strategy, the mean and the 90th 208 8 Figure 6: Distribution of city-wide average waiting time for treatment after the earthquake according to four emergency response plans, highlighting mean ( ) and 90th-percentile values ( P90%). The time is measured from when the patient is injured by the earthquake until he or she is treated in an operating room. Strategies 1 and 2 are baselines with limited coordination (LC) capacities, whereas strategies 3 and 4 introduce higher coordination (HC) capacities across the whole system level for resource allocation and patient transfers. The ambulance usage and the treatment spatial distribution are show in Figure 7 for each plan. percentile estimates of city-wide waiting time are 22 and 41 days, respectively, outperforming the 209 rst baseline strategy, but not signi cantly (Figure 6). With this strategy, the system mostly relies 210 on the largest two hospitals, located in the highlighted centric districts, to meet the demands of 211 surgical procedures. Figure 7b shows the corresponding distributions of treated patients and trans- 212 fers. The two largest hospitals treat 2.8k patients, 60% of the total demand, though their functional 213 operating rooms only constitute 44% of the total. Because multiple hospitals with non-functional 214 operating rooms send patients to the same large hospitals under this strategy, their operating rooms 215 over ow. Moreover, such an strategy leads to heavy use of roads from the periphery to the city 216 center. For example, our mean estimates indicate that 231 patients would have to be transported 217 from the southernmost hospital alone to the largest hospital, nearly twice as many as the maximum 218 number of transfers between any hospital pair in the rst baseline strategy. 219 1.5 Strategy 3: Sharing Ambulances 220 In the third strategy for e ective emergency response, hospitals transfer patients across the system 221 (see Methods). In addition, they share their ambulance resources across the system. This strategy 222 represents an emergency plan that requires high coordination at the system level. During the 223 emergency response following the next big earthquake in Lima, emergency managers could deploy 224 this policy as soon as they collect information on the actual status of the operating rooms in 225 hospitals and the actual distributions of injuries, which often takes a few days after the disaster 226 depending on its magnitude.40{43Using this strategy, the mean and 90-th percentile estimates of 227 city-wide waiting time are 10 and 19 days, respectively, signi cantly outperforming the rst and 228 second strategies by factors of 3 and 2.2 in the mean estimates, respectively (Figure 6). Because 229 9 (a) (b) (c) (d) Figure 7: Spatial distribution of patient treatment and transfers. These values represent the average transfers and treated people according to the uncertain distribution of earthquake casualties and residual capacities in the hospital system. (a) Strategy 1: hospitals with unavailable operating rooms can transfer patients to the closest working hospital. (b) Strategy 2: hospitals with unavailable operating rooms can transfer patients to the largest working hospital. (c) Strategy 3: ambulances are shared in the hospital system. (d) Strategy 4: EMTs deploy 15 additional mobile operating rooms during the emergency response. The gure only shows roads between hospital pairs that transferred at least 5 patients. 10 an optimization formulation is used at the system level under this policy, patient transfers are 230 e ective at transporting patients from lower-capacity zones to higher-capacity zones, leading to a 231 more e ective use of the functional operating rooms across the city. 232 The distribution of treated patients matches the distribution of the residual operating room 233 capacity in the system (Figure 7c). 2.6k people are treated in the hospitals in the four highlighted 234 centric districts, which represents 56% of the total patients in the city and closely approaches the 235 residual capacities in this zone, 59% of the total functional operating rooms. Unlike the second 236 strategy, this strategy does not overload the capacities in the two largest hospitals by sending 237 most patients to them, instead it distributes patients across the system according to the residual 238 capacities of each hospital. Additionally, unlike with the rst strategy, ambulance capacities do 239 not bottleneck the system with this policy. Because hospitals share all ambulances in the system, 240 the ambulances will work where they are most needed, from the periphery to the city center. An 241 emergency plan that implements such a strategy will lead to a more balanced use of ambulances 242 through the city and thus it ooads critical roads. With this strategy, close to 90 patients would 243 have to be transported from the southernmost hospital to the largest hospital, less than half the 244 number with the second baseline strategy. 245 1.6 Strategy 4: Deployment of Additional Operating Rooms by EMTs 246 In the fourth policy for e ective emergency response, EMTs will deploy 15 additional mobile operat- 247 ing rooms to alleviate high demand-capacity gaps across the system (see Methods). We assume that 248 the additional operating rooms will be functioning three days after the earthquake. This strategy 249 deploys the operating rooms in close proximity to existing hospitals to leverage their triage ar- 250 eas and additional resources such as personnel, power generators or backup water. As with the 251 third strategy, hospitals are also able to share ambulance capacities across the city. Because an 252 optimization formulation is also used, the additional operating rooms are e ectively deployed in 253 locations and quantities that are critical to improve the performance of the emergency response. 254 Using this policy, the mean and 90-th percentile estimates of city-wide waiting times are 8 and 15 255 days, respectively (Figure 6). As expected, these estimates outperform the response with the rst 256 policy due to the additional operating rooms in the system. 257 Our analysis strategically locates the additional operating rooms in the periphery, mainly in the 258 southernmost and northernmost zones (Figure 7d). By deploying eld hospitals in the periphery, 259 more patients can be treated there, ooading the hospitals in the center. With this policy, 2.3k 260 people, representing 48% of the total patients, are treated in the four previously highlighted centric 261 districts, 16% fewer patients than with the rst policy. As a result, fewer patients have to be 262 transported from the periphery to the city center, ooading critical roads even more. In this 263 case, the southernmost hospital will only have to transfer 76 patients to the largest hospital, a 17% 264 reduction than with the rst policy. If the EMTs deploy more eld hospitals using this methodology, 265 then treatment times will be further reduced, the periphery will be better supplied with needed 266 resources, and the usage of ambulances and critical roads will be further ooaded in the city. 267 2 Discussion 268 We present a methodology for characterizing the emergency response of hospital systems after 269 earthquakes and designing policies to treat patients e ectively. Our methodology establishes the 270 11 groundwork for assessing the value of hospital system coordination through a metric that measures 271 the performance of high-coordination emergency policies in terms of waiting times and e ective use 272 of hospital resources. Because our methodology considers hospitals in a large urban center behave 273 as a system, emergency managers and resilience ocers can apply our methodology to a whole 274 large city and evaluate optimal patient transfer strategies between hospitals, e ectively allocate 275 ambulances in the system, and guide the deployment of eld hospitals. 276 We found that a M 8.0 earthquake in Lima will cause a spatial distribution of casualties that 277 does not match the post-earthquake capacities of the hospital system (Figure 5). The zones with 278 higher post-earthquake capacity are located in the city center, in clear contrast with the zones 279 with higher post-earthquake demands of health services. Large numbers of patients are located in 280 the periphery, where, unlike the city center, de cient construction practices have rapidly increased 281 the seismic vulnerabilities of the housing infrastructure.44The neighborhoods in the periphery 282 tend to be populated by families with less income and wealth, so the disparities in disaster risk 283 overlap with the economic disparities.45This overlap will exacerbate the critical conditions in 284 the periphery because these families will have less resources to obtain treatment and medicine 285 from private hospitals, relying mostly on public hospitals. During an emergency response, this 286 uneven vulnerability pro le in the housing sector exacerbates the resource centralization problem 287 of hospitals in the city, thus leaving the neighborhoods in the periphery predominately underserved 288 during an emergency response. Because in many cities, neighborhoods in peripheral zones have 289 precarious access to health services46and high concentration of seismic vulnerabilities,47these 290 observations in Lima can be extrapolated to multiple urban centers in Latin America and even in 291 developed countries. 292 Emergency planners who aim to treat patients in the city e ectively must address these dis- 293 parities by either transporting patients from lower-capacity to higher-capacity zones or supplying 294 the lower-capacity zones with additional resources to meet demands of health services. Though 295 emergency managers can elaborate multiple reasonable strategies to implement such emergency re- 296 sponse measures, our ndings shows that strategies based on deeper coordination between hospitals 297 prove to be signi cantly more e ective than the ones with less ability to coordinate (Figure 6). 298 Turning to the two baseline strategies with reduced coordination, both of them have similar low 299 performance even though the system behavior di ers fundamentally. Because in the rst strategy 300 hospitals transfer patients to the next closest working hospital only when they do not have available 301 operating rooms, ambulance capacities are the bottleneck in some hospitals, whereas they are not 302 even used in other ones. Additionally, the rst strategy distributes the number of people treated 303 at each hospital roughly evenly (Figure 7a). Such a treatment distribution results in hospitals 304 with fewer resources treating similar patient numbers as hospitals with more resources, making the 305 system inecient and increasing waiting times. 306 The second strategy is also a baseline that enables hospitals to transfer patients only to the 307 hospital with the largest number of functional operating rooms in the city. Contrary to the rst 308 strategy, the second one distributes the number of treated patients highly unevenly (Figure 7b). 309 Because multiple hospitals often end up transferring patients to the two largest hospitals in the 310 city, these two hospitals largely over ow their capacities. In addition, the second strategy requires 311 heavy use of ambulance resources to transfer patients from the periphery to the city center, where 312 these two largest hospitals are located. 313 In contrast, the third and fourth strategies that have deeper coordination signi cantly improve 314 the emergency response performance because under these policies, hospitals share resources and 315 12 leverage strategic system-level information. Because ambulances are shared across the system, 316 they do not bottleneck the system as in rst baseline strategy and start to be e ectively used in 317 critical zones. Patients are also strategically transferred across the entire hospital system, leading 318 the spatial distribution of treated people matches the post-earthquake capacities of the hospital 319 system (Figure 7c). As a result, the citywide-average waiting times to treat patients decrease 320 quite signi cantly compared to the previous baseline strategies that only allowed reduced pairwise 321 coordination by factors larger than 2 (Figure 6). 322 In practice, emergency managers will need a few days after the next large earthquake to collect 323 necessary system-level information including both casualties and the residual functionality.40{43324 However, they can use the trends shown in Figure 7c to establish more informed policies for prag- 325 matic implementation. Emergency planners should use this type of assessment to create earthquake 326 preparedness plans. These plans may include implementing reliable communication lines between 327 hospitals more likely to transfer large number of patients. Additionally, preparedness plans may 328 include a usage prioritization of critical emergency corridors that are more likely to be heavily used 329 by ambulances or an identi cation of alternative roads in case these corridors have known seismic 330 vulnerabilities. 331 Another advantage of policies with deeper coordination the strategic deployment of EMTs and 332 their additional mobile operating rooms. Besides decreasing waiting times, EMTs that supply oper- 333 ating rooms to the most underserved zones also achieve a more ecient use of ambulances because 334 fewer transfers are needed. EMTs might nd it practical to locate these additional operating rooms 335 at the city center, where equipment mobilization is easier and availability of doctors and nurses 336 is higher. However, such a plan may lessen the ability of the hospital system to treat patients 337 because more patients would have to be transferred from the periphery to the center, overloading 338 the roads and potentially over owing ambulance capacities. Thus, robust earthquake preparedness 339 plans should be developed based on a thorough understanding of the uneven distribution of capac- 340 ity and demand of health services in an earthquake aftermath. E ective plans will capitalize on 341 the methodology and information provided here to better prepare cities facing high signi cant risk 342 from future large earthquakes. 343 Acknowledgements 344 We thank Abhinav Bindal, Jacqueline Li, and Dr. Maryia Markhvida for helping conceptualize the 345 rst version of the patient transfer analysis during an emergency response. We also thank Prof. 346 Sandra Santa Cruz, from the Ponti cia Universidad Cat olica del Per u for granting access to the 347 hospitals' seismic vulnerability information, and Prof. Carlos Zavala and Miguel Estrada from the 348 Centro Peruano-Japon es de Investigaciones S smicas y Mitigaci on de Desastres (CISMID) and the 349 Universidad Nacional de Ingenier a (UNI) for providing access to the seismic microzonation data in 350 Lima. Additionally, we thank Jill O'Nan for helpful revision of the paper writing. We acknowledge 351 the nancial support by the John A. Blume Fellowship from the Civil Engineering Department at 352 Stanford University. 353 13 3 Methods 354 3.1 Network Flow Model and Optimization Formulation 355 We model the post-earthquake hospital treatment process as a minimum cost time-varying net- 356 work ow (MCTVNF) problem.48,49In our MCTVNF formulation, a directed graph G= (N;E ) 357 represents the hospital system, where n=jNjis the number of graph nodes, and e=jEjis the 358 number of graph edges. We use a discrete time model with a nite time horizon tfwith time-steps 359 dt, thus the time t2T:f0;dt;2dt;:::;t fg. At each time t, each hospital has two nodes: one 360 triage node where patients are received into the hospitals, and one discharge node where patients 361 go after they complete their treatment. Each graph node is associated to an index iand a time 362 t, where hospitals' triage areas have indexes i2 :f1;2;:::;n hg, and the discharge areas have 363 indexesi2 :fnh+ 1;nh+ 2;:::; 2nhg, wherenhis the number of hospitals in the system. To 364 de ne a one-to-one correspondence within the indices of a hospital, if its triage index is i2, then 365 its discharge index is i+nh2. Figure 8 shows an example of a network representation at time 366 tfor a system with three hospitals, where the triage nodes are in red and the discharge nodes in 367 blue. In this model, the decision variables are both the ows through the edges and the patient 368 queues in thee triage nodes. These variables will track how many patients will stay in triage, be 369 treated or be transferred to other hospitals. 370 4 Dynamic Flipped b1(t) b2(t) b3(t)b4(t) b5(t) b6(t)x1,2(t)x2,1(t) x2,3(t) x3,2(t)x3,1(t)x3,1(t)x1,4(t) x2,5(t)x3,6(t) 3 Figure 8: System model with three hospitals at time tas a directed graph. The system model used for the application to Lima has 41 hospitals, i.e., 41 triage and 41 discharge nodes. Each graph node is associated to a time-variant demand-supply variable bi(t). In triage nodes, 371 bi(t) represents the number of people arriving to hospital, thus they are analyzed as source nodes 372 with nonnegative ows: bi(t)>= 0;8i2. In the discharge nodes, bi(t) represent the number 373 of patients who nish their treatment and exit the hospital at time t, thus they are analyzed as 374 sink nodes with nonpositive ows: bi(t)<= 0;8i2. We assume that patients who nish the 375 treatment process and exit the hospital do not to return to the hospital system during the time 376 horizontf. 377 Each graph edge is associated to a ow of patients xi;j(t) that leaves node iat timetto go to 378 14 nodej. In this formulation, edges fully connect the triage nodes to allow hospitals to redistribute 379 their patient loads to potentially any other hospital within the system according to their available 380 ambulances. Additionally, each triage node is connected to its respective discharge node to represent 381 the patient treatment process within a hospital. Figure 8 shows the edges and respective ows 382 between triage nodes from di erent hospitals and between triage and discharge nodes within same 383 hospitals for the system with three hospitals. At each time t, the owxi;j(t) has a maximum bound 384 ui;j(t) and a travel time i;j(t). In this discrete formulation, it is considered that the ow xi;j(t) 385 leaves the node iat timetand reaches the node jat timet+ i;j(t). 386 For the edges connecting triage areas, ui;j(t) represents the maximum number of patients who 387 can be transported from triage ito triagejin a di erent hospital according to the available 388 transportation resources (e.g., ambulances available in the hospital), and i;j(t) is the transportation 389 time of the patients from triage ito triagej. For the application to Lima, ui;j(t) and i;j(t) in 390 these edges were de ned according to the ambulance capacities in each hospital and the travel times 391 from pre-earthquake trac conditions, respectively. When vulnerability data for the transportation 392 system in Lima is available, our model will be able to leverage existing risk models for transportation 393 systems50,51to adjust travel times to post-earthquake trac conditions. 394 For the edges connecting triage nodes iwith their respective discharge nodes j=i+nh,ui;i+nh(t) 395 represents the maximum number of patients who can be treated according to the available medical 396 resources for the type and severity of the patients' injuries, and i;i+nh(t) is the treatment time. For 397 the application to Lima, ui;j(t) and i;j(t) in these edges were de ned according to the functional 398 operating rooms in each hospital and average treatment times in operating rooms in previous 399 earthquakes.52400 Additionally, we de ne yi(t) as a storage variable at each triage node to represent the patients 401 who wait in the hospital queue to either be treated within the hospital or be transported to another 402 hospital with more available resources. 403 3.1.1 Optimization of Performance Metric 404 We evaluate both waiting times and e ective use of ambulances as the system performance metric, 405 thus the metric includes two objective functions. The rst objective function C1(X) measures 406 waiting time across the city as the average time that a patient would take since the earthquake 407 until completing treatment in the operating room. 408 C1(X) =P t2TP i2; j=i+nhft+ i;jgxi;j(t)dt P t2TP i2bi(t)(1) Xrepresents a vector containing all the decision variables of ow xi;j(t) in edges and the storage 409 yi(t) in the triage nodes. The numerator of C1(X) represents the total number of patients passing 410 through the operating rooms (from each triage node i2 to the corresponding discharge node 411 j=i+nh2) multiplied by their respective times to complete treatment, whereas the denominator 412 is the total number of patients arriving to the triage areas. The time horizon tfis carefully chosen 413 to have enough modeling time to treat the patients. However, in few simulations with signi cant 414 number of patients and not many functional operating rooms, a couple of terms are added to the 415 numerator, one with the remainder patients in the triage areas, tfP i2yi(tf)dt, and another 416 with the remainder patients in the ambulances, tfP i2xi;i+nh(tf)dt, in order to properly 417 incorporate the unmet demands at the end of the simulation into C1(X). 418 15 The second objective function measures ambulance usage as the total number of patients trans- 419 ported in ambulances. The objective function C2(X) is normalized by the total number of patients 420 analogously to C1(X). 421 C2(X) =P t2TP i2P j2xi;j(t)dtP t2TP i2bi(t)(2) We de ne a system cost C(X) as a weighted sum of C1(X) andC2(X) to nd a Pareto-optimal 422 solution. 423 C(X) = 1C1(X) + 2C2(X) (3) After assessing multiple 1and 2values, we minimized C(X) using values of 0.90 and 0.1, 424 respectively. Smaller 2values resulted in inecient ambulance usage with small reductions in 425 waiting times, requiring some patients to be transferred multiple times in ambulances before being 426 treated. Larger 2signi cantly increased waiting times, thus these 2values do not appropriately 427 represent that the priority in the formulation is to minimize waiting times over to use ambulances 428 with eciency. We nd the best set of decisions ^X, vector that contains the values of ow variables 429 xi;j(t) and storage variables yi(t) which minimize C(X). 430 ^X= argminxi;j(t);yi(t)C(X) (4) The decision variables are subject to the constraints in Equations 5, 6, 7, and 8. Equation 5 431 represents patient ow conservation, which guarantees that all the patients coming into the hospital 432 system stay within the system until they leave through the discharge nodes. 433 xi;i+nh+X j2xi;j(t)X j2xj;i(t i;j(t)) +yi(t+dt)yi(t) =b(i);8i2;t2T (5) Equations 6 and 7 represent ow capacity constraints. Equation 6 ensures that the people in 434 the operating rooms do not exceed the unitary capacities ui;i+nh, whereui;i+nhis estimated as the 435 number of functional operating rooms in the hospital iover the number of surgeries per day. We 436 assumed that each surgery takes 4 hours, and that hospitals will be functional 24 hours during the 437 emergency response using multiple personnel shifts. Such treatment rate equals the rates in foreign 438 eld hospitals after the 2004 Indonesia earthquake/tsunami.52439 0xi;i+nh(t) ui;i+nh1;8i2;t2T (6) Equation 7 ensures that the patient transfers do not exceed the total unitary transportation 440 capacities in a hospital, where ui;jis the unitary capacity if all ambulances of a hospital were only 441 transferring patients from triage itoj.ui;jequals the number of ambulances in the hospital times 442 the number of patients transported per ambulance trip over the number of round trips that the an 443 ambulance can make from triage node itoj. We retrieved travel time information from Google 444 Maps API to estimate the round trip numbers and assumed that each ambulance trip can take up 445 to two patients. 446 0X j2xi;j(t) ui;j1;8i2;t2T (7) 16 Equation 8 ensures that the number of patients waiting in the hospitals' triage queues are 447 properly represented by a non-negative number. 448 0yi(t);8i2;t2T (8) Equations 6, 7 and 8 introduce a model relaxation. Whereas the number of patients who 449 are treated, transported or waiting in the queue can only be non-negative integers, the formulation 450 expands the variables' domain to include real numbers. This relaxation ensures that the formulation 451 is tractable. Thus, because the cost and the constraint functions are linear combinations of the 452 decision variables, we solve this minimization as a linear programming problem using the simplex 453 algorithm in GLPK of the cvxopt implementation in Python.53,54454 3.1.2 Model Adaptation for Baseline Strategies 1 and 2 455 Both baseline strategies have limited coordination capacity and only allow each hospital to transfer 456 patients to only one single hospital with functional operating rooms instead of multiple ones. Thus, 457 to represent these strategies, the model ignores multiple transfer edges in the ow model, reducing 458 the elements of the edge set E. In the rst baseline strategy, only the edges going from hospitals 459 without functional operating rooms to the closest hospitals are activated. In the second baseline 460 strategy, only the edges going from hospitals without functional operating rooms to the hospital 461 with the are largest number of functional operating rooms are activated. 462 Because the model is solved multiple times according to the number of patients and functional 463 operating rooms in the earthquake simulation, then the edge connectivity varies from simulation to 464 simulation. With strategies 1 and 2, the number of edges in the model is signi cantly reduced, thus 465 we modeled larger time horizons. We selected a time horizon Tfof 100 days, which is suciently 466 long period to treat treat all earthquake patients in most simulations, and a time step dtof 1 day. 467 3.1.3 Model Adaptation for Strategy 3: Sharing Ambulances 468 Strategy 3 does not need to disconnect edges in the model. Yet, it modi es the transportation 469 edges' capacity constraints to enable hospitals to share ambulance resources. Thus, the constraint 470 in Equation 7 is relaxed as follows. 471 0X i2aiX j2xi;j(t) pi;jX i2ai;8t2T (9) Equation 9 ensures that unitary transportation capacities are not exceeded at a system level at 472 each time step, where airepresents the number of ambulances of hospital i. All the other constraints 473 remain the same. Because modeling this policy requires higher edge connectivity than the baseline 474 strategies and thus has more computational demands, the time horizon tfwas reduced to 40 days. 475 It was veri ed that such a variation did not a ect the optimization because less modeling time was 476 needed as a result of shorter optimal waiting times with the strategies 3 and 4 (Figure 6). The 477 time stepdtwas kept equal to 1 day. 478 17 3.1.4 Model Adaptation for Strategy 4: Deployment of Additional Operating Rooms 479 by EMTs 480 Strategy 4 requires an additional modi cation to the constraint on the operating room capacity in 481 Equation 6. This strategy allows EMTS to increase hospital capacities by introducing additional 482 mobile operating rooms in close proximity to them as follows. 483 0xi;j(t)qi ui1;8i2;j=i+nh2;t2Tf0;dt;:::;t sg (10) Equation 10 ensures that hospitals can increase their unitary operation room capacities by qi 484 after the time tsat which the operating rooms in the eld hospitals are deployed in the city. In 485 addition the sum of the additional resources distributed across the system cannot exceed the total 486 capacityQsupplied by all the eld hospitals in the region as follows. 487 0X i2qiQ (11) All the other constraints remain the same. These modi cations barely change the optimization 488 complexity. Thus, we kept the time horizon equal to 40 days and the time step equal to 1 day. 489 3.2 Earthquake Casualty Modeling 490 We utilize an earthquake multiseverity casualty model previously developed by the authors15to 491 evaluate the spatial distribution of injuries requiring surgical treatment after the M 8.0 earthquake. 492 The model is probabilistic and uses ground shaking estimates to propagate the earthquake inten- 493 sity to building damage according to the building seismic vulnerability55and the site-speci c soil 494 conditions in Lima.56Next, the model uses information on building occupancy to provide proba- 495 bilistic estimates of the spatial distribution of injuries and fatalities in the city. The validity of the 496 model results was veri ed14by comparing the casualties and fatality levels in the city to empirical 497 formulas18and with fatality-to-collapse building data from the 2005 Pakistan earthquake.57498 The model categorizes injuries into three severities. The second- and third-degree severity re- 499 quire specialized medical attention and hospitalization, however, unlike the second degree, the third 500 one requires immediate rescue and treatment to avoid death.58{60We considered that 100% of the 501 patients with third-degree injuries, for example, having punctured organs or crush syndrome with 502 exposed wounds, plus 10% of patients with second-degree injuries, for example, having compound 503 bone fractures, will require surgical treatment in operating rooms. We considered that patients 504 arrive to the closest hospital during a period of 4 days after the earthquake in accordance to the 505 evidence from previous earthquakes.39,61Thus, in the ow model the demand-supply variable bi(t) 506 is larger than 0 in the triage nodes during the rst four days after the earthquake. We considered 507 that patients wait in triage zones to until an operating room is available in the hospital or until 508 they are transferred to other hospitals. 509 3.3 Seismic Analysis for hospital functionality 510 We utilize earthquake simulation to model the functionality of operating rooms during the emer- 511 gency response.62Hospitals are complex infrastructure, whose post-earthquake functionality de- 512 pends on multiple components: structural damage; damage in mechanical, electrical components 513 18 and medical equipment; utility failure; shortage of medical supplies (i.e., oxygen, blood), and short- 514 age of medical personnel.7,8,13,63Hospitals with slight structural damage can lose partial or total 515 functionality as a result of damage and loss of the other components of hospitals.64516 To capture these e ects, we analyzed that the structural vulnerability55of the +700 buildings 517 belonging to the 41 healthcare campuses in the city according to the earthquake shaking intensity 518 and the soil conditions on site. Then, we used a Bernoulli distribution to model loss of function- 519 ality that can occur due to failure of components di erent to the hospitals' structure according 520 the \Hospital Safety Index" (HSI). HSI is based on a qualitative evaluation of multiple hospital 521 components including buildings' nonstructural elements such as equipment and backup medical 522 resources, and technical and organization capacities in the hospitals' personnel.37HSI has three 523 categories: \A", \B" and \C", ranging from the best to the lowest performance. We used a di erent 524 Bernoulli distribution for each HSI category. We considered that operating rooms in buildings with 525 no structural damage have 1, 0.75, and 0.5 of functionality probability for categories \A", \B" and 526 \C", respectively, whereas that in buildings with slight structural damage, operating rooms have 527 0.6, 0.45, and 0.3 of functionality probability. Operating rooms in buildings with larger damage 528 levels were considered completely nonfunctional. 529 The 41 campuses in the dataset are part of the public healthcare system led by the Peruvian 530 Health Ministry (MINSA) and the Social Security (Essalud). Even though there is a growing private 531 healthcare system, most of the health care services are provided by the public system in Lima.65532 Physicians who work full time in the public healthcare system often work part-time in the private 533 system,66thus, in an emergency, they would aim to provide services in the public system rather 534 than in the private one. We consider that studying the response of the public sector represents a 535 robust starting point to characterize the earthquake emergency response of the hospital system in 536 Lima. 537 We supplemented the hospitals' building information with the number of ambulance in each 538 hospitals. Because, a few hospitals have no ambulances, we considered that during the emergency 539 response the local government or private institutions will supply one ambulance to each of these 540 hospitals so that each hospital is able to mobilize patients. 541 3.4 Earthquake Shaking 542 We studied the tectonics of the M 8.0 1940 earthquake and located the rupture area in the region 543 delimited by the earthquake aftershock zone.26We de ned the rupture dimensions along the fault 544 strike and dip directions using an empirical function based on subduction zone earthquake data.28545 Next, we evaluated the ground shaking in a grid of 1kmx1km using site-speci c lognormal dis- 546 triutions. We evaluated three ground shaking intensity measures, Peak ground acceleration, PGA , 547 spectral acceleration at 0.3s, Sa(0:3), and spectral acceleration at 1s, Sa(1:0s). We selected these 548 intensity measures to better capture the response of multiple typologies of buildings in the inven- 549 tory according to their predominant period of vibration. The log-mean and log-standard-deviation 550 values of the intensity measures were extracted from empirical formulas that relate magnitude, site 551 distance, and soil conditions to the ground shaking.29We included within-67and between-68event 552 correlations in the intensity measures. The between-event correlations introduce spatial correlations 553 to the ground shaking. 554 19 3.5 Data Availability 555 All the data to reproduce the ndings of the paper can be found at https://purl.stanford.edu/ 556 dp530wq8437 . 557 3.6 Code Availability 558 All the computer code to reproduce the ndings of the paper can also be found at https://purl. 559 stanford.edu/dp530wq8437 . 560 20 4 Supplementary Information 561 Lima is a fast-growing megacity with a population close to 10 million people.21Though the center 562 of the city is denser, the peripheral areas of the city have become heavily populated over the last 563 few decades. Currently, close to three million people live in peripheral zones in slums,69where 564 families are low-income, who often start constructing their homes with precarious materials, e.g., 565 wooden shacks, and then upgrade them to con ned-masonry buildings over timespans ranging from 566 a few years to decades.70Figure 9 shows how heavily populated the peripheries are. The popula- 567 tion distribution in this plot represents the average number of people over 24 hours in grids of 1 568 km2. Population density is dynamic, but often people spend most time at their residential build- 569 ings, mainly during nighttime. Thus, we considered that this average distribution is a reasonable 570 representation of nighttime population densities. 571 Figure 9: Spatial distribution of population density in Lima per km2. Data obtained from LandScan.71 The intervals in the two plots represent quintiles on the spatial data. The seismic analysis included the assessment of the the number of fatalities and injured people 572 with three types of severities caused by an M 8.0 earthquake occurring at nighttime in Lima, when 573 people are often within their houses. As observed in previous earthquakes, the model considers 574 that most casualties are caused by earthquake damage to buildings in the city. Our mean estimates 575 indicate that the M 8.0 earthquake will cause 60.3k people with injuries of severity 1, 18.9k of 576 severity 2, 2.8k of severity 3, and 5.6k immediate fatalities. People with injuries of severity 1 will 577 require basic medical aid and no hospitalization. People with injuries of severity 2 will require 578 hospital treatment, but the injuries are not life-threatening in the short term, and people with 579 severity 3 will require immediate hospitalization otherwise injuries become life-threatening.14,60580 21 Figures 10a and 10b show the mean spatial distributions of injured people with severity 2 and 3 in 581 the city. Because in the model, casualties are result of building damage, the spatial distribution of 582 patients with severity 2 and 3 are heavily cross-correlated and particularly concentrated in areas 583 with large number of buildings that collapse. 584 (a) (b) Figure 10: Spatial distribution of patients with a) injuries of severity 2 and b) injuries with severity 3. The intervals in the two plots represent quintiles on the spatial data. The seismic analysis also included the assessment of the ability of hospitals to function after the 585 earthquake. Figure 11 shows the mean functionality ratio of operating rooms in the 41 healthcare 586 campuses that were analyzed. The ratio represents how likely are operating rooms to function 587 after the M 8.0 earthquake considering their structural vulnerabilities, their \HSI" score, and their 588 proneness to experience large shaking intensities as a result of the proximity to the earthquake fault 589 or the soil conditions. Though the absolute number of functional operating rooms was higher at 590 the center of the city, the spatial patterns of of functionality ratios did show a strong prevalence 591 of high ratios in particular zones of the city. Instead of geographical location, construction year 592 was a better indicator of the hospitals' ability to function, as most newer hospitals showed higher 593 functionality ratios. 594 22 Figure 11: Spatial distribution post-earthquake functionality ratio. Newer hospitals, some of the located in the city center, tend to perform better than the older ones. Few hospitals that did not have operating room capacities under normal conditions were assigned 0% functionality ratio. 23