### Network flow model and optimization formulation

We model the post-earthquake hospital treatment process as a minimum cost time-varying network flow (MCTVNF) problem^{48,49}. In our MCTVNF formulation, a directed graph **G** = (**N**, **E**) represents the hospital system, where *n* = ∣**N**∣ is the number of graph nodes, and *e* = ∣**E**∣ is the number of graph edges. We use a discrete time model with a finite time horizon *t*_{f} with time-steps *d**t*, thus the time *t* ∈ **T**: {0, *d**t*, 2*d**t*, …, *t*_{f}}. At each time *t*, each hospital has two nodes: one triage node where patients are received into the hospitals, and one discharge node where patients go after they complete their treatment. Each graph node is associated to an index *i* and a time *t*, where hospitals’ triage areas have indexes *i* ∈ **Γ**: {1, 2, …, *n*_{h}}, and the discharge areas have indexes *i* ∈ **Λ**: {*n*_{h} + 1, *n*_{h} + 2, …, 2*n*_{h}}, where *n*_{h} is the number of hospitals in the system. To define a one-to-one correspondence within the indices of a hospital, if its triage index is *i* ∈ **Γ**, then its discharge index is *i* + *n*_{h} ∈ **Λ**. Fig. 7 shows an example of a network representation at time *t* for a system with three hospitals, where the triage nodes are in red and the discharge nodes in blue. In this model, the decision variables are both the flows through the edges and the patient queues in thee triage nodes. These variables will track how many patients will stay in triage, be treated or be transferred to other hospitals.

Each graph node is associated to a time-variant demand-supply variable *b*_{i}(*t*). In triage nodes, *b*_{i}(*t*) represents the number of people arriving to hospital, thus they are analyzed as source nodes with non-negative flows: *b*_{i}(*t*) > = 0, ∀ *i* ∈ **Γ**. In the discharge nodes, *b*_{i}(*t*) represent the number of patients who finish their treatment and exit the hospital at time *t*, thus they are analyzed as sink nodes with non-positive flows: *b*_{i}(*t*) < = 0, ∀ *i* ∈ **Λ**. We assume that patients who finish the treatment process and exit the hospital do not to return to the hospital system during the time horizon *t*_{f}. Each graph edge is associated to a flow of patients *x*_{i,j}(*t*) that leaves node *i* at time *t* to go to node *j*. In this formulation, edges fully connect the triage nodes to allow hospitals to redistribute their patient loads to potentially any other hospital within the system according to their available ambulances. In addition, each triage node is connected to its respective discharge node to represent the patient treatment process within a hospital. Fig. 7 shows the edges and respective flows between triage nodes from different hospitals and between triage and discharge nodes within same hospitals for the system with three hospitals. At each time *t*, the flow *x*_{i,j}(*t*) has a maximum bound *u*_{i,j}(*t*) and a travel time *τ*_{i,j}(*t*). In this discrete formulation, it is considered that the flow *x*_{i,j}(*t*) leaves the node *i* at time *t* and reaches the node *j* at time *t* + *τ*_{i,j}(*t*).

For the edges connecting triage areas, *u*_{i,j}(*t*) represents the maximum number of patients who can be transported from triage *i* to triage *j* in a different hospital according to the available transportation resources (e.g., ambulances available in the hospital), and *τ*_{i,j}(*t*) is the transportation time of the patients from triage *i* to triage *j*. For the application to Lima, *u*_{i,j}(*t*) and *τ*_{i,j}(*t*) in these edges were defined according to the ambulance capacities in each hospital and the travel times from pre-earthquake traffic conditions, respectively. When vulnerability data for the transportation system in Lima is available, our model will be able to leverage existing risk models for transportation systems^{50} to adjust travel times to post-earthquake traffic conditions.

For the edges connecting triage nodes *i* with their respective discharge nodes *j* = *i* + *n*_{h}, ({u}_{i,i+{n}_{h}}(t)) represents the maximum number of patients who can be treated according to the available medical resources for the type and severity of the patients’ injuries, and ({tau }_{i,i+{n}_{h}}(t)) is the treatment time. For the application to Lima, *u*_{i,j}(*t*) and *τ*_{i,j}(*t*) in these edges were defined according to the functional operating rooms in each hospital and average treatment times in operating rooms in previous earthquakes^{51}.

In addition, we define *y*_{i}(*t*) as a storage variable at each triage node to represent the patients who wait in the hospital queue to either be treated within the hospital or be transported to another hospital with more available resources.

### Optimization of performance metric

We evaluate both waiting times and effective use of ambulances as the system performance metric, thus the metric includes two objective functions. The first objective function *C*_{1}(*X*) measures waiting time across the city as the average time that a patient would take since the earthquake until completing treatment in the operating room.

$${C}_{1}({bf{X}})=frac{{sum }_{tin {bf{T}}}{sum }_{begin{array}{c}iin {boldsymbol{Gamma }},\ j = i+{n}_{{rm{h}}}end{array}}{t+{tau }_{i,j}}times {x}_{i,j}(t)times dt}{{sum }_{tin {bf{T}}}{sum }_{iin {boldsymbol{Gamma }}}{b}_{i}(t)}$$

(1)

**X** represents a vector containing all the decision variables of flow *x*_{i,j}(*t*) in edges and the storage *y*_{i}(*t*) in the triage nodes. The numerator of *C*_{1}(**X**) represents the total number of patients passing through the operating rooms (from each triage node, *i* ∈ **Γ** to the corresponding discharge node *j* = *i* + *n*_{h} ∈ **Λ**) multiplied by their respective times to complete treatment, whereas the denominator is the total number of patients arriving to the triage areas. The time horizon *t*_{f} is carefully chosen to have enough modeling time to treat the patients. However, in few simulations with significant number of patients and not many functional operating rooms, a couple of terms are added to the numerator, one with the remainder patients in the triage areas, *t*_{f} × ∑_{i∈Γ}*y*_{i}(*t*_{f}) × *d**t*, and another with the remainder patients in the ambulances, ({t}_{{rm{f}}}times {sum }_{iin {boldsymbol{Gamma }}}{x}_{i,i+{n}_{{rm{h}}}}({t}_{{rm{f}}})times dt), in order to properly incorporate the unmet demands at the end of the simulation into *C*_{1}(**X**).

The second objective function measures ambulance usage as the total number of patients transported in ambulances. The objective function *C*_{2}(**X**) is normalized by the total number of patients analogously to *C*_{1}(**X**).

$${C}_{2}({bf{X}})=frac{{sum }_{tin {bf{T}}}{sum }_{iin {boldsymbol{Gamma }}}{sum }_{jin {boldsymbol{Gamma }}}{x}_{i,j}(t)times dt}{{sum }_{tin {bf{T}}}{sum }_{iin {boldsymbol{Gamma }}}{b}_{i}(t)}$$

(2)

We define a system cost *C*(**X**) as a weighted sum of *C*_{1}(**X**) and *C*_{2}(**X**) to find a Pareto-optimal solution.

$$C({bf{X}})={alpha }_{1}times {C}_{1}({bf{X}})+{alpha }_{2}times {C}_{2}({bf{X}})$$

(3)

After assessing multiple *α*_{1} and *α*_{2} values, we minimized *C*(**X**) using values of 0.90 and 0.1, respectively. Smaller *α*_{2} values resulted in inefficient ambulance usage with small reductions in waiting times, requiring some patients to be transferred multiple times in ambulances before being treated. Larger *α*_{2} significantly increased waiting times, thus these *α*_{2} values do not appropriately represent that the priority in the formulation is to minimize waiting times over to use ambulances with efficiency. We find the best set of decisions (hat{{bf{X}}}), vector that contains the values of flow variables *x*_{i,j}(*t*) and storage variables *y*_{i}(*t*) which minimize *C*(**X**).

$$hat{{bf{X}}}={{rm{argmin}}}_{{x}_{i,j}(t);{y}_{i}(t)}quad C({bf{X}})$$

(4)

The decision variables are subject to the constraints in Equations (5), (6), (7), and (8). Equation (5) represents patient flow conservation, which guarantees that all the patients coming into the hospital system stay within the system until they leave through the discharge nodes.

$${x}_{i,i+{n}_{{rm{h}}}}+{sum }_{jin {boldsymbol{Gamma }}}{x}_{i,j}(t)-{sum }_{jin {boldsymbol{Gamma }}}{x}_{j,i}(t-{tau }_{i,j}(t))+{y}_{i}(t+dt)-{y}_{i}(t)=b(i),quad forall iin {boldsymbol{Gamma }},tin {bf{T}}$$

(5)

Equations (6) and (7) represent flow capacity constraints. Equation (6) ensures that the people in the operating rooms do not exceed the unitary capacities ({u}_{i,i+{n}_{{rm{h}}}}), where ({u}_{i,i+{n}_{{rm{h}}}}) is estimated as the number of functional operating rooms in the hospital *i* over the number of surgeries per day. We assumed that each surgery takes 4 hours, and that hospitals will be functional 24 hours during the emergency response using multiple personnel shifts. Such treatment rate equals the rates in foreign field hospitals after the 2004 Indonesia earthquake/tsunami^{51}.

$$0le frac{{x}_{i,i+{n}_{{rm{h}}}}(t)}{{u}_{i,i+{n}_{{rm{h}}}}}le 1,quad forall iin {boldsymbol{Gamma }},tin {bf{T}}$$

(6)

Equation (7) ensures that the patient transfers do not exceed the total unitary transportation capacities in a hospital, where *u*_{i,j} is the unitary capacity if all ambulances of a hospital were only transferring patients from triage *i* to *j*. *u*_{i,j} equals the number of ambulances in the hospital times the number of patients transported per ambulance trip over the number of round trips that an ambulance can make from triage node *i* to *j*. We retrieved travel time information from Google Maps API to estimate the round trip numbers and assumed that each ambulance trip can take up to two patients.

$$0le {sum }_{jin Gamma }frac{{x}_{i,j}(t)}{{u}_{i,j}}le 1,quad forall iin {boldsymbol{Gamma }},tin {bf{T}}$$

(7)

Equation (8) ensures that the number of patients waiting in the hospitals’ triage queues are properly represented by a non-negative number.

$$0le {y}_{i}(t),quad forall iin {boldsymbol{Gamma }},tin {bf{T}}$$

(8)

Equations (6), (7), and (8) introduce a model relaxation. Whereas the number of patients who are treated, transported or waiting in the queue can only be non-negative integers, the formulation expands the variables’ domain to include real numbers. This relaxation ensures that the formulation is tractable. Thus, because the cost and the constraint functions are linear combinations of the decision variables, we solve this minimization as a linear programming problem using the simplex algorithm in GLPK of the cvxopt implementation in Python^{52}.

### Model adaptation for baseline strategies 1 and 2

Both baseline strategies have limited coordination capacity and only allow each hospital to transfer patients to only one single hospital with functional operating rooms instead of multiple ones. Thus, to represent these strategies, the model ignores multiple transfer edges in the flow model, reducing the elements of the edge set **E**. In the first baseline strategy, only the edges going from hospitals without functional operating rooms to the closest hospitals are activated. In the second baseline strategy, only the edges going from hospitals without functional operating rooms to the hospital with the largest number of functional operating rooms are activated.

Because the model is solved multiple times according to the number of patients and functional operating rooms in the earthquake simulation, then the edge connectivity varies from simulation to simulation. With strategies 1 and 2, the number of edges in the model is significantly reduced, thus we modeled larger time horizons. We selected a time horizon *t*_{f} of 100 days, which is sufficiently long period to treat all earthquake patients in most simulations, and a time step *d**t* of 1 day.

### Model adaptation for sharing ambulances

Strategy 3 does not need to disconnect edges in the model. Yet, it modifies the transportation edges’ capacity constraints to enable hospitals to share ambulance resources. Thus, the constraint in Equation (7) is relaxed as follows.

$$0le sum _{iin {boldsymbol{Gamma }}}{a}_{i}sum_{jin {boldsymbol{Gamma }}}frac{{x}_{i,j}(t)}{{p}_{i,j}}le sum_{iin {boldsymbol{Gamma }}}{a}_{i},quad forall tin {bf{T}}$$

(9)

Equation (9) ensures that unitary transportation capacities are not exceeded at a system level at each time step, where *a*_{i} represents the number of ambulances of hospital *i*. All the other constraints remain the same. Because modeling this policy requires higher edge connectivity than the baseline strategies and thus has more computational demands, the time horizon *t*_{f} was reduced to 40 days. It was verified that such a variation did not affect the optimization because less modeling time was needed as a result of shorter optimal waiting times with the strategies 3 and 4 (Fig. 5). The time step *d**t* was kept equal to 1 day.

### Model adaptation for deployment of more operating rooms

Strategy 4 requires an additional modification to the constraint on the operating room capacity in Equation (6). This strategy allows EMTS to increase hospital capacities by introducing additional mobile operating rooms in close proximity to them as follows.

$$0le frac{{x}_{i,j}(t)-{q}_{i}}{{u}_{i}}le 1,quad forall iin {boldsymbol{Gamma }},j=i+{n}_{h}in {boldsymbol{Lambda }},tin {bf{T}}-{0,dt,ldots ,{t}_{{rm{s}}}}$$

(10)

Equation (10) ensures that hospitals can increase their unitary operation room capacities by *q*_{i} after the time *t*_{s} at which the operating rooms in the field hospitals are deployed in the city. In addition, the sum of the additional resources distributed across the system cannot exceed the total capacity *Q* supplied by all the field hospitals in the region as follows.

$$0le sum_{iin {boldsymbol{Gamma }}}{q}_{i}le Q$$

(11)

All the other constraints remain the same. These modifications barely change the optimization complexity. Thus, we kept the time horizon equal to 40 days and the time step equal to 1 day.

### Earthquake casualty modeling

We utilize an earthquake multiseverity casualty model previously developed by the authors^{17} to evaluate the spatial distribution of injuries requiring surgical treatment after the M 8.0 earthquake. The model is probabilistic and uses ground shaking estimates to propagate the earthquake intensity to building damage according to the building seismic vulnerability^{53} and the site-specific soil conditions in Lima^{54}. Next, the model uses information on building occupancy to provide probabilistic estimates of the spatial distribution of injuries and fatalities in the city. The validity of the model results was verified^{16} by comparing the casualties and fatality levels in the city to empirical formulas^{55} and with fatality-to-collapse building data from the 2005 Pakistan earthquake^{56}.

The model categorizes injuries into three severities. The second- and third-degree severity require specialized medical attention and hospitalization, however, unlike the second degree, the third one requires immediate rescue and treatment to avoid death^{57}. We considered that 100% of the patients with third-degree injuries, for example, having punctured organs or crush syndrome with exposed wounds, plus 10% of patients with second-degree injuries, for example, having compound bone fractures, will require surgical treatment in operating rooms. We considered that patients arrive to the closest hospital during a period of 4 days after the earthquake in accordance to the evidence from previous earthquakes^{32}. Thus, in the flow model the demand-supply variable *b*_{i}(*t*) is larger than 0 in the triage nodes during the first four days after the earthquake. We considered that patients wait in triage zones to until an operating room is available in the hospital or until they are transferred to other hospitals.

### Seismic analysis for hospital functionality

We utilize earthquake simulation to model the functionality of operating rooms during the emergency response^{58}. Hospitals are complex infrastructure, whose post-earthquake functionality depends on multiple components: structural damage; damage in mechanical, electrical components and medical equipment; utility failure; shortage of medical supplies (i.e., oxygen, blood), and shortage of medical personnel^{9,10,14,59}. Hospitals with slight structural damage can lose partial or total functionality as a result of damage and loss of the other components of hospitals^{60}.

To capture these effects, we analyzed that the structural vulnerability^{53} of the +700 buildings belonging to the 41 healthcare campuses in the city according to the earthquake shaking intensity and the soil conditions on site. Then, we used a Bernoulli distribution to model loss of functionality that can occur owing to failure of components different to the hospitals’ structure according the Hospital Safety Index (HSI). HSI is based on a qualitative evaluation of multiple hospital components including buildings’ nonstructural elements such as equipment and backup medical resources, and technical and organization capacities in the hospitals’ personnel^{29}. HSI has three categories: A, B, and C, ranging from the best to the lowest performance. We used a different Bernoulli distribution for each HSI category. We considered that operating rooms in buildings with no structural damage have 1, 0.75, and 0.5 of functionality probability for categories A, B, and C, respectively, whereas that in buildings with slight structural damage, operating rooms have 0.6, 0.45, and 0.3 of functionality probability. Operating rooms in buildings with larger damage levels were considered completely non-functional.

The 41 campuses in the data set are part of the public healthcare system led by the Peruvian Health Ministry (MINSA) and the Social Security (Essalud). Even though there is a growing private healthcare system, most of the healthcare services are provided by the public system in Lima^{61}. Physicians who work full time in the public healthcare system often work part-time in the private system^{62}, thus in an emergency, they would aim to provide services in the public system rather than in the private one. We consider that studying the response of the public sector represents a robust starting point to characterize the earthquake emergency response of the hospital system in Lima.

We supplemented the hospitals’ building information with the number of ambulance in each hospitals. Because, a few hospitals have no ambulances, we considered that during the emergency response the local government or private institutions will supply one ambulance to each of these hospitals so that each hospital is able to mobilize patients.

### Earthquake shaking

We studied the tectonics of the M 8.0 1940 earthquake and located the rupture area in the region delimited by the earthquake aftershock zone^{23}. We defined the rupture dimensions along the fault strike and dip directions using an empirical function based on subduction zone earthquake data^{63}.

Next, we evaluated the ground shaking in a grid of 1 km × 1 km using site-specific lognormal distributions. We evaluated three ground shaking intensity measures, peak ground acceleration spectral acceleration at 0.3s, *S**a*(0.3*s*), and spectral acceleration at 1s, *S**a*(1.0*s*). We selected these intensity measures to better capture the response of multiple typologies of buildings in the inventory according to their predominant period of vibration. The log-mean and log-standard-deviation values of the intensity measures were extracted from empirical formulas that relate magnitude, site distance, and soil conditions to the ground shaking^{64}. We included within-^{65} and between-^{66}event correlations in the intensity measures. The between-event correlations introduce spatial correlations to the ground shaking.

### Model limitations and future work

Our formulation and case study advance the field of earthquake emergency response. However, there are existing limitations that must be addressed in future research. First, we model patient transfers that heavily rely on the proper functioning of the transportation system during the emergency response. However, our paper does not capture that the transportation network can be disrupted owing to potential bridge failures, seismic liquefaction, or infrastructure debris caused by the earthquake. Although there are models to capture these disruptions in a large city, we did not capture them because they require extensive seismic vulnerability data currently not available for Lima. Only a few large urban centers worldwide have sufficient vulnerability data to start capturing disruptions in the entire transportation network, e.g., the Bay Area in California^{50}. Because collecting such data often requires significant resources from cities, strategic surveying of critical elements (e.g., critical bridges and roads) in the system can be an effective starting point for cities that lack sufficient vulnerability data. Our results can be used to support such a strategic surveying approach in large urban centers Instead of exhaustively collecting the vulnerability data of all components of the transportation network, the de facto approach in risk analysis, researchers can use our results to prioritize surveying the critical roads for patient transfers to evaluate its potential disruptions owing to structural failure, soil liquefaction, and debris blockage and propose strategic risk mitigation measures for the network. Similarly, the power, water, and sewage networks^{67} are pivotal for the proper functioning of the hospital facilitates. Our results can help future research focus on the power substations, power distribution lines, and water and sewage pipes serving the critical hospitals in the system. Focusing on these critical components will enable large urban centers which lack data to have a robust starting point for including the interactions between the hospital system and other critical urban systems.

Second, our case study only included an evaluation of an earthquake occurring at nighttime because we lacked data with the spatial distribution of the urban population at other times of the day. This spatial distribution is an important factor in the assessment of the spatial distributions of earthquake casualties because it enables researchers to track what vulnerable buildings are occupied. In our nighttime scenario, people are mostly at residential buildings, which are particularly vulnerable in Lima. We used LandScan information, which is available worldwide^{68}, as a proxy for the nighttime distribution of the population^{16,17}. Earthquake scenarios occurring during the daytime will have a different distribution of population and therefore a different spatial distribution of earthquake casualties. Although we were not able to track the time variations of urban population in Lima, existing work has already shown that mobile network data can be effectively used for it, e.g., Beijing in China^{69}. As more mobility data are rapidly collected by public and private institutions, plenty of opportunities to create large databases with time-variant high-resolution urban densities will be available in the near future. Deployments of models to track these time-variant urban densities will enable researchers to leverage the formulation presented here and extend the emergency response application to different times of the day in multiple cities so that emergency managers are better informed on how to protect their communities.

### Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.