Throughout history, pandemics have significantly influenced societal structures, economic stability, political systems, and impacted mortality rates.Reference Conti 17 In recent decades, the ongoing threat of emerging and reemerging pandemics and epidemics has posed a significant challenge to humanity.Reference Brodeur, Gray and Islam 11 In a similar yet distinct manner, the overall number of wars and conflicts worldwide has shown a concerning increase recently, following a relatively peaceful period after World War II.Reference Cunen, Hjort and Nygard 18 , Reference Davies, Pettersson and Oberg 19
Pandemics and military operations have long been interlinked due to the sociological and operational dynamics associated with armies and warfare.Reference Erkoreka 22 , Reference Flecknoe, Wakefield and Simmons 30 The levels of hygiene, self-care, and immune system stress in individuals provided fertile conditions for the outbreak of diseases within armies and the communities interacting with them.Reference Aligne 5 Moreover, the density of military bases and army deployments facilitates the spread of pathogens within the population.Reference Aligne 5 Notable examples include the 1917-1918 Spanish Influenza, which affected populations worldwide and particularly impacted the US military during World War I, and the COVID-19 outbreak in Ukraine during the Russia-Ukraine war starting in 2022.Reference Anghel and Jones 6 , Reference Erkoreka 23
The intersection of pandemics and military operations has placed economic and health care systems worldwide under considerable strain.Reference Krueger, Uhlig and Xie 44 , 80 , Reference Wu, Perrings and Kinzig 93 To prepare for such events, researchers frequently design policies and plans for a wide range of scenarios using mathematical models.Reference Ferguson, Laydon and Nedjati-Gilani 26 , Reference Fernandes and Maldonado 27 , Reference Lazebnik and Blumrosen 48 , Reference Lemaitre, Grantz and Kaminsky 54 , Reference Eylul Oruc, Baxter and Keskinocak 64 Models related to warfare are often less accessible to the general public and, consequently, can be difficult to locate in the literature.Reference Sengar and Lohiya 76 Unlike, numerous models have been proposed to simulate the spread of pandemics.Reference Lazebnik, Bunimovich-Mendrazitsky and Shaikhet 50 , Reference Rahimi, Chen and Gandomi 71 , Reference Vytla, Ramakuri and Peddi 87 For example,Reference Miller 60 proposed an SIR (Susceptible-Infected-Recovered) based disease spread model that incorporates multiple transmission pathways, including direct contact and vectors. This model effectively captures the dynamics of pandemic spread across large populations,Reference Chen 14 and introduced an extended SIR model that differentiates between symptomatic and asymptomatic cases of COVID-19, assigning distinct infection paths for each. In this model, susceptible individuals can become asymptomatic immediately upon infection, while those who develop symptoms first enter an exposed stage. Yan et al. (2018) proposed an extended spatio-temporal SIR model for the spread of sexually transmitted diseases on a bipartite random contact network.Reference Yan, Zhang and Ma 94
Building on top of these modeling efforts, a growing body of work has focused on the capacity of health care systems to support the treatment of infected individuals during large-scale pandemics.Reference Cascini, Hoxhaj and Zaçe 13 For instance, Kozyreff et al. (2021) utilized an SIR model incorporating the dynamics of hospitalized patients, assumed to be proportional to the number of infected individuals, with considerations for delay and exponential decay during their stay in health care centers.Reference Kozyreff 42 Notably, the hospitalized population does not affect other epidemiological states, which is crucial for accurately modeling health care center dynamics. This model was applied to COVID-19 pandemic data from the first 3 months in 5 different regions: Belgium, France, Italy, Switzerland, and New York City (US). Similarly, Huang et al. (2021) introduced a risk-stratified SIR-HCD model that incorporates variables representing the dynamics of low- and high-contact subpopulations, with shared parameters to manage their respective rates over time.Reference Huang, Chu and Shams 35 The authors tested their model using data on daily reported hospitalizations and cumulative mortality from COVID-19 in Harris County, TX.
Despite the extensive epidemiological-mathematical research, a model that captures the dynamics of a pandemic during wartime and its influence on the health care system’s ability to provide necessary services has yet to be developed. This gap leaves policymakers without crucial insights for planning in such scenarios. Hence, in this work, we propose a novel spatio-temporal extended SIR-based model that includes civilian and military health care systems with asymmetric administration policies to study the effect of pandemics on warfare performance from the health care perspective. In addition, we provide an agent-based simulation implementation of the proposed model with a reinforcement learning model that is able to approximate optimal patient administration policies for the entire health care system, including both war and civilian zones to reduce overall mortality (from both war and pandemic), operating as an applicative tool to investigate efficient administration polices with different properties without the need to manually search for them. A schematic view of the study’s structure is presented in Figure 1.
A schematic view of the study’s design.
First, we formalize the model by introducing 4 key components: spatial movement dynamics, pandemic spread dynamics, hospitalization dynamics, and warfare dynamics. Next, we utilize agent-based simulation to solve the model and optimize patient administration policies using reinforcement learning. Finally, we establish a realistic setup to evaluate the pandemic- and war-related deaths, analyzing the impact of various patient administration strategies on these outcomes.
The remaining paper is organized as follows: Section 2 presents the related work for pandemic spread during wars and spatio-temporal modeling practices for pandemic spread. Section 3 describes the proposed model’s mathematical formalization. Section 4 outlines the implementation of the proposed model as a computer simulator and optimal patient administration using a reinforcement learning model. Section 5 provides a comprehensive evaluation of the proposed model. Finally, Section 6 provides a conclusion on the possible applications and limitations followed by suggestions for future work.
Related Work
Warfare has been studied from various perspectives, including sociological,Reference Kuznar 45 political,Reference Balcells and Steele 8 and economic,Reference Liadze, Macchiarelli and Mortimer-Lee 56 among others. In addition, over the last century, advancements in biology, medicine, and engineering have greatly enhanced our epidemiological understanding.Reference Piret and Boivin 68 As such, the relationship between the 2 started to emerge in the literature and a growing body of work captured the interplay between warfare and pandemics. In this section, we review the unique relationship between warfare and pandemic spread. Furthermore, we present recent modeling practices for pandemic spread which will be the mathematical base for our model.
Pandemic During War
In the context of large-scale pandemics and wars, one can detect 2 central and relatively recent global-scale events - the Spanish influenza pandemic of 1917-1918 during the First World WarReference Erkoreka 23 and the COVID-19 pandemic during the Russian-Ukrainian war.Reference Anghel and Jones 6
Following the sequence of events, we first examine the Spanish influenza pandemic, which stands out as one of the deadliest pandemics in recent history.Reference Edwards 21 Estimates place the global death toll between 21.5 million and 39.3 million, with some scholars suggesting even higher numbers, though these may be exaggerated.Reference Patterson and Pyle 67 The initial wave of the pandemic, occurring in the spring of 1918, was relatively mild with few fatalities. However, a drastic change occurred in the summer of the same year when the virus underwent a mutation, becoming highly virulent and leading to millions of deaths worldwide in October and November.Reference Ansart, Pelat and Boelle 7 A less severe third wave followed in early 1919, while the fourth and final wave extended into the initial months of 1920. The demographic most affected comprised young, healthy adults aged 15-44. Mortality rates varied across countries and continents, with estimates in Europe ranging from 1.1%-1.2%.Reference Erkoreka 24 Importantly, Oxford et al. proposed that the pandemic originated at the British military base in Étaples, located in northern France in the Pas-de-Calais department.Reference Oxford, Lambkin and Sefton 65 Significant during World War I, this base housed 100 000 soldiers within a 12 square kilometer area. Adjacent to the base were coastal marshes frequented by migratory birds. In the vicinity, numerous farms kept pigs, ducks, and geese for the soldiers’ consumption, along with horses used for transportation. The convergence of dense soldier populations, various livestock, and the extensive use of 24 types of war gases, several of which were mutagenic, could potentially explain the initial outbreak of the epidemic between December 1916 and March 1917. Moreover, Erkoreka et al. suggested that the Spanish influenza was exported by boat to the US from Tanger and the Spanish colonies in the north of Africa, highlighting the influence of military personnel movement as a vector of pandemic spread.Reference Erkoreka 23
The Russian-Ukraine war began in February 2022, during the third year of the COVID-19 pandemic. By that time, the pandemic, caused by the SARS-CoV-2 pathogen, had resulted in over 6 million deaths.Reference Choudhary, Saied and Priyanka 15 Before the pandemic, Ukraine’s health care system had many shortcomings, making it difficult to detect and control the spread of COVID-19.Reference Choudhary, Saied and Priyanka 15 These challenges included a lack of essential health care equipment, inadequate personal protective gear, a shortage of skilled health care workers, limited availability of infectious disease departments across different regions of Ukraine, insufficient surveillance measures, and a scarcity of authorized testing laboratories. These issues indicate that Ukraine was already struggling with the COVID-19 crisis, a situation made worse by the ongoing conflict.Reference Choudhary, Saied and Priyanka 15 Due to the Russian invasion of Ukraine, approximately 5 million Ukrainians immigrated to neighboring European countries, such as Poland, Romania, Russia, Hungary, Moldova, Slovakia, and Belarus. The invasion also caused population mixing inside of the country as people relocated from war zones.Reference Chumachenko and Chumachenko 16 During the invasion, the Ukrainian community experienced major disruptions in primary health care services and faced challenges accessing basic preventive medicine.Reference Patel, Moncayo and Conroy 66 Additionally, COVID-19 vaccination efforts and routine immunizations nearly ceased nationwide. The rise in COVID-19 cases in Ukraine caused concern for the EU. The large influx of refugees—the largest since World War II—posed a risk of spreading the virus beyond Ukraine’s borders, potentially triggering new outbreaks across Europe.Reference Quinn, Dhabalia and Roslycky 70
Spatio-temporal Pandemic Modeling
Pandemics can arise from various sources, including sexually transmitted diseases,Reference Quinn 69 socially influenced behavior,Reference Djillali, Bentout and Touaoula 20 and airborne diseases.Reference Ferguson, Cummings and Fraser 25 Airborne pandemics have raised the most concern due to their high infection rates and the wide population they can affect without restriction. This category of airborne pandemics includes various pathogens such as influenza, Lassa virus, COVID-19, and Nipah virus, among others.Reference Wang, Prather and Sznitman 88
Modeling pandemic spread and control is a complex interdisciplinary endeavor that requires understanding various domains like mathematics, biology, medicine, sociology, and economics; all of which influence or are influenced by pandemics and pandemic intervention policies (PIPs).Reference Lazebnik, Bunimovich-Mendrazitsky and Ashkenazi
49
, Reference Rock, Brand and Moir
72 These models can generally be categorized into statistical models and mechanistic models. Statistical models rely on data-driven approaches without specific assumptions on dynamics, often using statistical and machine learning methods for forecasting some properties related to the pandemic.Reference Ivorra, Ferrandez and Vela-Perez
36
, Reference Salgotra, Gandomi and Gandomi
74 On the other hand, mechanistic models, commonly exemplified by the Susceptible-Infected-Recovered (SIR) model, operate based on theoretical principles to explain pandemic dynamics.Reference Fernández-Villaverde and Jones
28 For instance, the SIR model assumes that the population is divided into susceptible, infected, and recovered individuals such that each infected individual infected on average
$ \beta $
susceptible individuals and that infected individuals recover and become recovered at a rate
$ \gamma $
.Reference Kermack and McKendrick
39 The SIR model’s simplicity often falls short in capturing real-world pandemic dynamics,Reference Moein, Nickaeen and Roointan
61 leading to the development of numerous extensions and adaptations.Reference Adiga, Dubhashi and Lewis
1
, Reference Wiratsudakul, Suparit and Modchang
91 These extensions are typically divided into spatial and temporal categories. Spatial extensions are represented either by norm-based or graph-based models, addressing the well-mixed population assumption used in the SIR model. Temporal extensions introduce more types of interactions and properties for individuals in the population, which play a role in the spread of the pandemic.Reference Lazebnik
46
Recently, several spatio-temporal extended SIR-based models have been proposed to better capture the spread of the COVID-19 pandemic.Reference Bearman, Moody and Stovel 9 , Reference Keeling 37 , Reference Keeling and Eames 38 , Reference Lazebnik and Alexi 47 , Reference Lazebnik, Shami and Bunimovich-Mendrazitsky 52 For instance, Viguerie et al. proposed a partial differential equation (PDE) with the SEIRD (E - exposed, D - dead) model based spatio-temporal.Reference Viguerie, Lorenzo and Auricchio 86 Different diffusion parameters are also factored in for various population groups with spatial movement across a large population modeled using an inhomogeneous random walk, converging to a second-order differential operator in the limit. The authors applied this model specifically to the COVID-19 outbreak in Lombardy, Italy, focusing on the initial 6 months of the pandemic. Goel and Sharma (2021) introduced a location graph-based SIR model where each graph node represents a location inhabited by a subset of the population.Reference Goel and Sharma 32 In their model, susceptible individuals can become infected either through contact with infectious individuals at the same node (local infection) or through infectious individuals from neighboring nodes who visit their node (global infection). The authors also factored in social connectivity, which represents the average number of individuals moving between nodes, thus facilitating global infection. In our model, we follow this trend, proposing a graph-based extended SIR model for the pandemic spread.
Moreover, multiple models also take into consideration the complex relationship between health care systems and pandemics.Reference Hespanha, Chinchilla and Costa 34 , Reference Weissman, Crane-Droesch and Chivers 90 For example, Fox et al. (2022) used a stochastic extended SIRD model with symptomatic infectious, asymptomatic infectious, and hospitalized individuals such that the model’s parameters are binomial distributed rather than constant.Reference Fox, Lachmann and Tec 31 The authors used data from COVID-19 in Texas (US) to predict the number of hospitalized individuals over time due to the pandemic. Similarly, Li and Zhang (2017) integrated the average quality of treatment into the SIR model to capture the pandemic spread with different levels of health care systems.Reference Li and Zhang 55 In our model, we also include a health care system, extending it to handle both pandemic-related patients and soldiers with war injuries.
Model Definition
During times of war, countries are divided into civilian and war zones,Reference Verwimp, Justino and Bruck 85 each hosting distinct yet interconnected populations of civilians and soldiers. In civilian zones, individuals engage in various activities such as commuting for work, accessing social services, or leisure purposes like tourism.Reference Findley and Young 29 Similarly, war zones consist of multiple locations where soldiers are deployed and often commute between, either for operational duties or to civilian areas for rest and recuperation. The effectiveness of soldiers in fulfilling their duties in war zones is critically dependent on adequate personnel availability.Reference Simonton 78 Injured soldiers in war zones are evacuated to health care centers—either military or civilian—where the proximity of the facility to the site of injury significantly impacts the quality of care and recovery outcomes.Reference McGuirk and Porell 58 , Reference Nicholl, West and Goodacre 62 Compounding the challenges of war, pandemics can emerge, affecting both soldiers and civilians. Infected individuals may present as asymptomatic or symptomatic, with the latter at greater risk of death without timely medical intervention.Reference He, Guo and Mao 33 To improve survival rates, symptomatic individuals often require hospitalization, emphasizing the importance of accessible health care during such crises.Reference Lechien, Chiesa-Estomba and Place 53 , Reference Wu, Li and Shi 92 Figure 2 presents a schematic view of the pandemic spread with a dual health care system during a war model.
A schematic view of the pandemic spread with a dual health care system during war model. The dynamics can be divided into 4 interconnected components: pandemic spread (
$ {SEI}^a{I}^s RD $
model), war-related death and wounded soldiers (Lanchester model), commute between locations (graph-based model), and health care system (functional model).
To capture the epidemiological dynamics during warfare involving civilian and military populations within a dual-use health care system, we propose a spatio-temporal framework that integrates an extended SIR model for epidemiological dynamics and the Lanchester model for warfare dynamics. For simplicity, we divide the proposed model into 4 components: spatial movement, pandemic spread, hospitalization, and warfare. Based on these components, we assemble a single mathematical framework.
Formally, we define the spatio-temporal model that captures these dynamics as follows. A model,
$ M $
, is defined by the tuple
$ M\hskip0.2em := \hskip0.2em \left(G,P,H,W\right) $
such that
$ G $
is the graph of locations and possible commute between them,
$ P $
is the pandemic spread dynamics,
$ H $
is the association of health care centers with locations and their treatment dynamics, and
$ W $
is the war-related dynamics. Below, we define each component independently and then combine them into a single framework.
Spatial graph representation
Let us consider a graph
$ G\hskip0.2em := \hskip0.2em \left({V}_w,{V}_c,E\right) $
where
$ {V}_w $
is the set of war zone locations,
$ {V}_c $
is the set of civilian locations, and
$ E $
is the set of possible commutes between all locations, such that
$ E\hskip0.3em \subseteq \hskip0.3em \left({V}_w\cup {V}_c\right)\times \left({V}_w\cup {V}_c\right)\times {\mathrm{\mathbb{R}}}^2 $
. Each location contains 2 sub-populations - a civilian and a military. In terms of the pandemic spread, the populations are well-mixedFootnote
1 and heterogeneous. For the commute dynamics, each edge
$ e\hskip0.3em := \hskip0.3em \left(i,j.{c}_w,{c}_c\right)\in E $
between nodes
$ i $
and
$ j $
, has 2 weights
$ {c}_w $
and
$ {c}_c $
indicating the average commute between the 2 locations for military and civilian sub-populations, respectively. In practice,
$ {c}_w $
and
$ {c}_c $
are time-dependent, however, for simplicity, we use the average value over time.
Pandemic spread dynamics
We adopted the
$ {SEI}_a{I}_s RD $
model.Reference Cai and Zhou
12
, Reference Lazebnik, Bunimovich-Mendrazitsky and Shami
51
, Reference Nurlaila, Hidayat and Pardamean
63 Namely, the model considers a constant population with a fixed number of individuals
$ N $
. For simplicity, and given the short time horizon of interest, we abstract from population growth. Each individual belongs to 1 of the 6 groups: susceptible
$ \Big(S $
), exposed
$ E $
, asymptomatic infected
$ \Big({I}^a $
), symptomatic infected
$ \Big({I}^s $
), recovered
$ (R) $
, and dead
$ (D) $
such that
$ N=S+E+{I}^a+{I}^s+R+D $
. Individuals in the first group have no immunity and are susceptible to infection. When an individual in the susceptible group (
$ S $
) is exposed to the pathogen, the individual is transferred to the exposed state (
$ E $
) at rate
$ \beta $
. The individual stays in the exposed state for
$ \psi $
time steps in time and then transforms to either the asymptomatic infected group (
$ {I}^a $
) or symptomatic infected group (
$ {I}^s $
) with rate
$ \rho $
and
$ 1-\rho $
, respectively. The individual stays in the symptomatic infected group on average
$ {\gamma}_s $
time steps, after which the individual is transferred to the recovered group (
$ R $
) or the dead group
$ (D) $
with rate
$ \lambda $
and
$ 1-\lambda $
, respectively. All asymptomatic infected individuals stay in the infected
$ {I}^a $
group on average
$ {\gamma}_a $
time steps, after which the individual is transferred to the recovered group
$ (R) $
. The recovered are again healthy, no longer contagious, and immune from future infection. The epidemiological dynamics are formally described below using a system of ordinary differential equations:
$$ {\displaystyle \begin{array}{c}\frac{dS(t)}{dt}=-\beta S(t)\left({I}_s(t)+{I}_a(t)\right),\\ {}\frac{dE(t)}{dt}=\beta S(t)\left({I}_s(t)+{I}_a(t)\right)-\psi E(t),\\ {}\frac{dI_s(t)}{dt}=\rho \psi E(t)-{\gamma}_s{I}_s(t),\\ {}\frac{dI_a(t)}{dt}=\left(1-\rho \right)\psi E(t)-{\gamma}_a{I}_a(t),\\ {}\frac{dR(t)}{dt}={\gamma}_a{I}_a(t)+{\lambda \gamma}_s{I}_s(t),\\ {}\frac{dD(t)}{dt}=\left(1-\lambda \right){\gamma}_s{I}_s(t).\end{array}} $$
Figure 3 presents a schematic view of the pandemic spread models, divided into its epidemiological states.
A schematic view of the pandemic spread model presenting the epidemiological states and the transformations between them.
Hospitalization dynamics
The health care system is based on a set of health care service provider locations. Each health care service provider’s location is able to treat a number of individuals (
$ {\mu}_s $
) while providing the best clinical treatment possible which reduces the probability of dying due to the pandemic or war injuries (
$ \lambda, \varepsilon $
) in a rate
$ {\zeta}_s $
. However, once the number of individuals administrated to the health care center extends
$ {\mu}_s $
, the treatment performance is linearly decreasing until
$ {\mu}_f $
individuals are administrated to the health care center where the treatment performance is reaching
$ {\zeta}_f $
. In other words,
$ \Psi :\mathrm{\mathbb{N}}\to \left[0,1\right] $
, indicates the reduction in the probability of dying due to the pandemic (
$ \lambda $
) as a function of the number of administrated individuals:
$$ \Psi (x):= \left\{\begin{array}{l}{\zeta}_s\;\mathrm{if}\;x<{\mu}_s\\ {}{\zeta}_s-\frac{\zeta_f={\zeta}_s}{\mu_f-{\mu}_s}\left(x-{\zeta}_s\right)\;\mathrm{if}\;{\mu}_s\le x\le {\mu}_f\\ {}0\;\mathrm{if}\;x>{\mu}_f\end{array}\right. $$
Figure 4 presents a schematic view of the hospitalization performance as a function of the number of administrated individuals.
A schematic view of the hospitalization performance as a function of the number of administrated individuals.
War dynamics
To capture the war dynamics, we adopted the Lanchester modelReference Kress
43 and added wounded soldiers dynamics. Formally, we assume that there are only 2 combating sides and that the 2 sides comprise a homogeneous force. Let us consider
$ A(t) $
and
$ B(t) $
to be the sizes of the 2 armies at time
$ t $
. Each side has an average performance of reducing the size of the other army denoted by
$ {\varepsilon}_a $
and
$ {\varepsilon}_b $
, respectively. From the population that was removed from each army at some point in time
$ t $
, a portion is only wounded with injury rates of
$ {\theta}_a $
and
$ {\theta}_b $
, respectively. Wounded soldiers do not participate in the war and die out at a rate
$ d $
if not treated. Formally, this dynamic takes the form:
$$ {\displaystyle \begin{array}{c}\frac{dA(t)}{dt}=-{\varepsilon}_bB(t),\\ {}\frac{dB(t)}{dt}=-{\varepsilon}_aA(t),\\ {}\frac{dW^a(t)}{dt}={\theta}_b{\varepsilon}_bB(t)-{dW}^a(t),\\ {}\frac{dW^b(t)}{dt}={\theta}_a{\varepsilon}_aA(t)-{dW}^b(t).\end{array}} $$
Assembling the components into a single framework
To assemble the 4 components into a single framework, one has to outline the relationship between the components. First, the graph-based spatial component introduces a spatial separation to local interaction - either for the war or pandemic dynamics as well as a commute between the locations, it enforces the later dynamics to occur on each node of the graph and includes movement of the population. Next, the pandemic’s spread should consider both soldiers and civilians to fit the war dynamics and the commute dynamics that occur in parallel. The pandemic and war occurring for both armies in each location are combined into a single process. Finally, both symptomatic infected individuals and wounded soldiers are administrated to health care service providers according to some policy.
To this end, the framework is presented at the node level such that the dynamics in all the nodes occur in parallel over time. Because there are both civilian and war-related nodes, we describe the dynamics in each separately. Starting with the civilian node, the dynamics take the following form:
$$ {\displaystyle \begin{array}{c}\frac{dS^c(t)}{dt}=-\beta {S}^c(t)\left({I}_s^c(t)+{I}_a^c(t)+{I}_s^m(t)+{I}_a^m(t)\right)+ Commute\left({S}^c\right),\\ {}\hskip-8em \frac{dE^c(t)}{dt}=\beta {S}^c(t)\left({I}_s^c(t)+{I}_a^c(t)+{I}_s^m(t)+{I}_a^m(t)\right)-\psi {E}^c(t)\\ {}\hskip-22em + Commute\left({E}^c\right),\\ {}\hskip-2em \frac{dI_s^c(t)}{dt}=\rho \psi {E}^c(t)-{\gamma}_s{I}_s^c(t)+ Commute\left({I}_s^c\right)- hosp\left({H}^c(t),{I}_s^c(t)\right),\\ {}\hskip-12em \frac{dI_a^c(t)}{dt}=\left(1-\rho \right)\psi {E}^c(t)-{\gamma}_a{I}_a^c(t)+ Commute\left({I}_a^c\right),\\ {}\hskip-5em \frac{dR^c(t)}{dt}={\gamma}_a{I}_a^c(t)+{\lambda \gamma}_s{I}_s^c(t)+\left(\lambda +\Psi \left({H}^c(t)+{H}^m(t)\right)\right){H}^c(t)\\ {}\hskip-21em + Commute\left({R}^c\right),\\ {}\hskip-5em \frac{dD^c(t)}{dt}=\left(1-\lambda \right){\gamma}_s{I}_s^c(t)+\left(1-\lambda -\Psi \left({H}^c(t)+{H}^m(t)\right)\right){H}^c(t),\\ {}\hskip-18em \frac{dH^c(t)}{dt}= hosp\left({H}^c(t),{I}_s^c(t)\right)-{\gamma}_s{H}^c(t),\\ {}\frac{dS^m(t)}{dt}=-\beta {S}^m(t)\left({I}_s^c(t)+{I}_a^c(t)+{I}_s^m(t)+{I}_a^m(t)\right)+ Commute\left({S}^m\right),\\ {}\hskip-5em \frac{dE^m(t)}{dt}=\beta {S}^m(t)\left({I}_s^c(t)+{I}_a^c(t)+{I}_s^m(t)+{I}_a^m(t)\right)-\psi {E}^m(t)\\ {}\hskip-19em + Commute\left({E}^m\right),\\ {}\frac{dI_s^m(t)}{dt}=\rho \psi {E}^m(t)-{\gamma}_s{I}_s^m(t)+ Commute\left({I}_s^m\right)- hosp\left({H}^m(t),{I}_s^m(t)\right),\\ {}\hskip-10em \frac{dI_a^m(t)}{dt}=\left(1-\rho \right)\psi {E}^m(t)-{\gamma}_a{I}_a^m(t)+ Commute\left({I}_a^m\right),\\ {}\hskip-5em \frac{dR^m(t)}{dt}={\gamma}_a{I}_a^m(t)+{\lambda \gamma}_s{I}_s^m(t)+\left(\lambda +\Psi \left({H}^c(t)+{H}^m(t)\right)\right){H}^m(t)\\ {}\hskip-21.5em + Commute\left({R}^m\right),\\ {}\hskip-18em \frac{dH^m(t)}{dt}= hosp\left({H}^m(t),{I}_s^m(t)\right)-{\gamma}_s{H}^m(t),\end{array}} $$
where
$ n\hskip0.3em := \hskip0.3em \mid {V}_w\mid +\mid {V}_c\mid $
is the total number of nodes,
$ {m}_{i,j} $
is the average commute between nodes
$ i $
and
$ j $
,
$ Commute(X)\hskip0.2em := {\sum}_{j=1,j\ne k}^n\left({m}_{j,k}{X}_j(t)\right)-{\sum}_{j=1,j\ne k}^n\left({m}_{k,j}X(t)\right) $
, and
$ hosp\left(x,y\right) $
is the patient administration policy for symptomatically infected individuals.
For the war nodes, the dynamics occur for 2 army groups as it is assumed no citizens of either side of the war are present in such locations. Thus, the dynamics take the form:
$$ {\displaystyle \begin{array}{c}\frac{dS^a(t)}{dt}=-{\beta}^a{S}^a(t)\left({I}_s^a(t)+{I}_a^a(t)\right)-{\beta}^{ab}\left({I}_s^b(t)+{I}_a^b(t)\right)-{\varepsilon}_b\frac{S^a(t){\Lambda}_b(t)}{\Lambda_a(t)}+ Commute\left({S}^a\right),\\ {}\frac{dE^a(t)}{dt}={\beta}^a{S}^a(t)\left({I}_s^a(t)+{I}_a^a(t)\right)+{\beta}^{ab}\left({I}_s^b(t)+{I}_a^b(t)\right)-\psi {E}^a(t)-{\varepsilon}_b\frac{E^a(t){\Lambda}_b(t)}{\Lambda_a(t)}+ Commute\left({E}^a\right),\\ {}\frac{dI_s^a(t)}{dt}=\rho \psi {E}^a(t)-{\gamma}_s{I}_s^a(t)+ Commute\left({I}_s^a\right)- hosp\left({H}^a(t),{I}_s^a(t),{W}^a(t)\right),\\ {}\frac{dI_a^a(t)}{dt}=\left(1-\rho \right)\psi {E}^a(t)-{\gamma}_a{I}_a^a(t)-{\varepsilon}_b\frac{I_a^a(t){\Lambda}_b(t)}{\Lambda_a(t)}+ Commute\left({I}_a^a\right),\\ {}\frac{dR^a(t)}{dt}={\gamma}_a{I}_a^a(t)+{\lambda \gamma}_s{I}_s^a(t)-{\varepsilon}_b\frac{R^a(t){\Lambda}_b(t)}{\Lambda_a(t)}+\left(\lambda +\Psi \left({H}^a(t)\right)\right){\gamma}_s{H}^a(t)+\Psi \left({H}^a(t)\right){\gamma}_w{H}^a(t)+ Commute\left({R}^a\right),\\ {}\frac{dW^a(t)}{dt}={\theta}_b{\varepsilon}_b{\Lambda}_b(t)-{dW}^a(t)- hosp\left({H}^a(t),{I}_s^a(t),{W}^a(t)\right),\\ {}\frac{dD^a(t)}{dt}=\left(1-\lambda \right){\gamma}_s{I}_s^a(t)-\Big(1-\lambda -\Psi \left({H}^a(t)\right){\gamma}_s{H}^a(t),\\ {}\frac{dH^a(t)}{dt}= hosp\left({H}^a(t),{I}_s^a(t),{W}^a(t)\right)-{\gamma}_s{H}^a(t)-{\gamma}_w{H}^a(t),\\ {}\frac{dS^b(t)}{dt}=-{\beta}^b{S}^b(t)\left({I}_s^b(t)+{I}_a^b(t)\right)-{\beta}^{ab}\left({I}_s^a(t)+{I}_a^a(t)\right)-{\varepsilon}_a\frac{S^b(t){\Lambda}_a(t)}{\Lambda_b(t)}+ Commute\left({S}^b\right),\\ {}\frac{dE^b(t)}{dt}={\beta}^b{S}^b(t)\left({I}_s^b(t)+{I}_a^b(t)\right)+{\beta}^{ab}\left({I}_s^b(t)+{I}_a^b(t)\right)-\psi {E}^b(t)-{\varepsilon}_a\frac{E^b(t){\Lambda}_a(t)}{\Lambda_b(t)}+ Commute\left({E}^b\right),\\ {}\frac{dI_s^b(t)}{dt}=\rho \psi {E}^b(t)-{\gamma}_s{I}_s^b(t)+ Commute\left({I}_s^b\right),\\ {}\frac{dI_a^b(t)}{dt}=\left(1-\rho \right)\psi {E}^b(t)-{\gamma}_a{I}_a^b(t)-{\varepsilon}_a\frac{I_a^b(t){\Lambda}_a(t)}{\Lambda_b(t)}+ Commute\left({I}_a^b\right),\\ {}\\ {}\\ {}\\ {}\\ {}\hskip-15.3em \frac{dR^b(t)}{dt}={\gamma}_a{I}_a^b(t)+{\lambda \gamma}_s{I}_s^b(t)-{\varepsilon}_a\frac{R^b(t){\Lambda}_a(t)}{\Lambda_b(t)}+\left(\lambda +\Psi \left({H}^b(t)\right)\right){\gamma}_s{H}^b(t)+\Psi \left({H}^b(t)\right){\gamma}_w{H}^b(t)+ Commute\left({R}^b\right),\\ {}\frac{dW^b(t)}{dt}={\theta}_a{\varepsilon}_a{\Lambda}_a(t)-{dW}^b(t)- hosp\left({H}^b(t),{I}_s^b(t),{W}^b(t)\right),\\ {}\frac{dD^b(t)}{dt}=\left(1-\lambda \right){\gamma}_s{I}_s^b(t)-\Big(1-\lambda -\Psi \left({H}^b(t)\right){\gamma}_s{H}^b(t),\\ {}\frac{dH^b(t)}{dt}= hosp\left({H}^b(t),{I}_s^b(t),{W}^b(t)\right)-{\gamma}_s{H}^b(t)-{\gamma}_w{H}^b(t),\end{array}} $$
where
$ {\beta}^a $
,
$ {\beta}^b $
, and
$ {\beta}^{ab} $
are the infection rates of the first army, second army, and between armies. We denote
$ {\Lambda}_a(t)={S}^a(t)+{E}^a(t)+{I}_a^a(t)+{R}^a(t) $
and
$ {\Lambda}_b(t)={S}^b(t)+{E}^b(t)+{I}_a^b(t)+{R}^b(t) $
.
Computer Simulation
In this section, we present an agent-based simulation (ABS) approach to solve the proposed model (Eq. [4-5]), followed by a reinforcement learning-based model to derive a (near-)optimal patient administration policy.
System Solver Using Agent-based Simulation
We implemented the proposed model (Eq. [4-5]) using the ABS approach following the scheme proposed by Labeznik (2023).Reference Lazebnik
46 Formally, let us assume a population of agents allocated to a graph of location (
$ G $
) in some pre-defined distribution. These agents move and interact in discrete finite time steps
$ t\in \left[1,\dots, T\right] $
, where
$ T<\infty $
. In order to use the ABS approach, one has to define the agents in the dynamics as well as their 3 types of interactions: agent-agent, agent-environment, and spontaneous (i.e., depends only on the agent’s state and time).Reference Tesfatsion
82 To this end, for our model, each agent is represented by a timed finite state machineReference Alagar and Periyasamy
3 with its epidemiological-clinical status (
$ S,E,{I}^c,{I}^a,R,D,W,H $
), current location (
$ v\in {V}_m\cup {V}_c $
), and sociological status (solder or civilian). The simulation starts by generating the location graph and allocating the agents inside it with their walking dynamics, as indicated by the average commute rate between each of the 2 locations, and epidemiological state. Next, iteratively, until a stop condition is met or a pre-defined number of steps in time is reached, the simulation solves the dynamics of each node in a random order. First, the portion of the population in each node moves to other nodes according to the average commute rate between locations. Afterward, the epidemiological and war-related dynamics take place (Eq. [4-5]). Finally, individuals are allocated to health care centers according to some patient administration policies. Notably, health care centers and locations can be associated in many-to-many manner.
Patient Administration Policy Optimization
In this section, we outline a deep reinforcement learning (DRL) model to determine the optimal patient administration policy aimed at minimizing overall mortality, considering deaths from both war-related causes and the pandemic. The RL method leverages the agent-based simulation (ABS) model described earlier to iteratively improve decision-making regarding the allocation of patients to health care centers.
The goal of the DRL method is to identify a patient administration policy
$ \pi $
that minimizes the cumulative death count over the simulation period
$ t\in \left[0,T\right] $
. The administration policy
$ \left(\pi \right) $
determines the actions (patient administration) based on the current state of the system, which includes the epidemiological-clinical status of individuals, their locations, and the capacity and status of health care centers. Formally, the RL model’s state space at time
$ t $
(
$ {s}_t $
) is a comprehensive representation of the current situation, including, the epidemiological status of individuals (
$ S,E,{I}^s,{I}^a,R,D,W,H $
), their current location (
$ v\in {V}_m\cup {V}_c $
), sociological status, and current capacity and occupancy of health care centers, as well as the state of the other army in the war zones. The action space,
$ A $
, is composed of the set of patient administrations from all nodes in the graph to all health care centers and represented by a tensor of size
$ \mid {V}_m\cup {V}_c\mid \times \mid C\mid \times 2 $
such that
$ \mid C\mid $
indicates the number of health care centers where each value in the tensor is a 2-dimensional vector indicating the number of civilians and soldiers allocated to each of the health care centers. The reward function of the RL model for time
$ t $
is
$$ R\left({s}_t,{a}_t\right)=-\sum \limits_{v\in {V}_m\cup {V}_c}\left({D}_v(t)+{W}_v(t)\right). $$
We employ a deep Q-network (DQN) algorithm [92], a model-free, online, and off-policy RL method, to optimize the patient administration policy. A DQN algorithm is a value-based RL method that trains a critic to estimate the expected discounted cumulative long-term reward when following the optimal administration policy. For the neural network architecture, we used 3 hidden layers including an LSTM with 64 neurons and 2 fully connected layer of sizes 64 and 32 with ReLu activations. The output layer’s size corresponds to the number of health care centers (
$ \mid C\mid $
). Notably, during the training phase, actions involving nodes that cannot allocate civilians and/or soldiers to specific hospitals, or to hospitals that are already at full capacity, were masked by assigning a large negative value to their corresponding entries in the action space (
$ -{10}^9 $
).
The RL algorithm is implemented within the ABS framework, enabling the simulation of various patient administration strategies and their impacts on overall mortality. Through iterative learning, the administration policy
$ \pi $
converges towards an optimal strategy that effectively allocates patients to minimize deaths from both war-related causes and the pandemic. Formally, in order to allow the RL model to learn a representative policy, we repeat each ABS configuration multiple times and then alter the ABS configuration. Table 1 presents the parameter values of the DRL model. The parameter values are chosen using a manual trial-and-error to achieve a good administration policy’s performance.
Hyperparameters for training the Q-learning model
Evaluation
To evaluate the proposed model, we first established realistic, while synthetic, scenarios. The sociological and movement parameters such as the population size and number of nodes in the graphs (which represent central cities) aim to capture a large-scale region up to the entire western country.Reference Weber and Sciubba 89 For the pandemic’s spread, we focused on the COVID-19 pathogen due to the large amount of accurate recorded data associated with this pathogen.Reference Betthäuser, Bach-Mortensen and Engzell 10 The warfare and hospitalization parameters ranged greatly in the literature. As such, we adopted a range without anomalies. Table 2 summarizes the parameters and their values as these are used by the proposed model.
The model’s parameters’ description, value ranges, and sources
Initially, we train the proposed DRL model to learn a (near-)optimal policy for patient administration. Figure 5 shows the relative performance of the DRL model, with aims to decrease the total number of deaths from both the pandemic and war (see Eq. [6]), with respect to the normalized random allocation (baseline) that is computed for
$ n=100 $
random simulations. Initially, the DRL’s performance is comparable to the baseline one, as it is also initialized with random logic. Next, for the first 80 simulations, the DRL performance was worse than the baseline, up to 35% in comparison. This initial decline in performance can be attributed to the exploration phase of the DRL model which searches for policies that may be worse than random allocation, such as no-allocation at all. However, after this point, a relatively smooth exponential improvement in the DRL’s performance occurs, converging to around 30% of the baseline performance (i.e., 70% improvement) after 800 simulations. Notably, around the 300-350 and 625-675 simulations, a local divergence takes place which can be associated with a local optimum of the optimization process converging and requires more exploration to overcome it and re-converge to more global optima.Reference Zhang, Koppel and Zhu
97
The DRL’s administration policy’s average performance normalized to the naive administration policy for
$ n=100 $
simulations as a function of the number of training simulations.
Next, to evaluate the performance of the DRL model, we compared it to 2 baseline policies and 2 single-usage health care polices. Figure 6 shows the overall death as the mean
$ \pm $
standard deviation of
$ {n}_{test}=100 $
simulations following a
$ {n}_{train}=1000 $
simulation to train the DRL model for each case, divided into 4 patient administration policies - no patient administration, random allocation up to full capacity, optimized based on Eq. (6), DRL model optimized to reduce the death of soldiers alone, and DRL model optimized to reduce the death of civilians alone. The optimal configuration is found to be statistically significantly better based on ANOVA (Analysis of Variance) with Tukey post-hoc testReference St and Wold
79 with
$ p<0.05 $
. Unsuspectingly, the case with no health care system is significantly worse than the other cases, with slightly over 4% death of the population. The naive random allocation results in less than half of the percent of death while also the standard deviation is smaller, as indicated by the error bars. Moreover, comparing the single-usage health care cases, in the form of the soldiers and civilians alone cases, the latter is statistically significantly worse, with
$ p<0.01 $
based on a Mann–Whitney U test.Reference McKnight and Najab
59
Overall death as a function of different administration policies. The results are shown as the mean
$ \pm $
standard deviation of
$ n=100 $
simulations.
Figure 6 revealed that providing health care services to soldiers, on average, results in overall less death in a population. Nevertheless, taking civilians into account provides a statistically significant improvement. Thus, we further explore the death and injured dynamics of soldiers over time. Figure 7 shows the wounded soldiers and overall deaths (both soldiers and civilians) over time under the optimal DRL administration policy. The results are presented as the mean
$ \pm $
standard deviation of
$ n=100 $
simulations, assuming a fixed population size of
$ N=5\cdot {10}^6 $
, with 10% of the population comprising army personnel. For the injured soldiers, an initial linear increase is observed up to the 35th simulation step, aligning with the expectations of the Lanchester model (see Eq. [3]). Afterward, a more chaotic pattern emerges, characterized by a slower rate of increase, which can be attributed to treated soldiers returning from the health care system. This behavior persists until approximately the 75th simulation step. These dynamics are followed by another increase up to the 110th simulation step, after which the values stabilize for approximately 35 simulation steps. Subsequently, a sharp decline is observed, reaching zero around the 185th simulation step. This decline can be explained by the end of the war as the other army is eliminated. In a complementary manner, the total deaths exhibit a monotonically increasing trend with a sigmoid-like pattern, reflecting the combined fatalities of injured soldiers who succumbed to war-related injuries and both soldiers and civilians who died due to the pandemic. The impact of injured soldiers on total deaths is twofold: their presence drives a faster-than-linear increase in deaths, delayed relative to their numbers. This strain on the health care system, occupied with treating injured soldiers, limits its capacity to manage infected individuals, exacerbating pandemic spread and mortality.
The wounded soldiers and overall deaths (both soldiers and civilians) over time using the optimal DRL administration policy. The results are shown as the mean
$ \pm $
standard deviation of
$ n=100 $
simulations assuming a fixed population size of
$ N=5\cdot {10}^6 $
, with 10% of the population are in the army.
Figure 7 reveals non-linear and phase-dependent dynamics, which further increase the uncertainty of the dynamics. Because the model assumes know-in-advance war results, in the form of the Lanchester model, the uncertainty can be associated with the pandemic spread, which “breaks” the symmetry in the Lanchester dynamics at the node level. Thus, the event optimization horizon for which the DRL model is trained should result in different performances. To this end, Figure 8 shows the mean percent improvement over the random administration (baseline) policy for
$ {n}_{test}=100 $
simulations followed by
$ {n}_{train}=1000 $
simulations. Notably, the model’s performance improves when the event optimization horizon aligns with the average recovery rate. Interestingly, this trend holds even for a 28-day event optimization horizon, suggesting that other dynamics, such as warfare, have a secondary effect on administration policy optimization due to their greater predictability.
The mean percent improvement over the random administration (baseline) policy for
$ {n}_{test}=100 $
simulations followed by
$ {n}_{train}=1000 $
simulations.
Conclusion
In this study, we have developed a novel spatio-temporal model to address the complex dynamics of dual-use health care administration policy during a pandemic in wartime conditions. By extending the SIR model to include the effects of both civilian and military populations, we offer a comprehensive framework that captures the intricate interactions between pandemic spread, war dynamics, and health care system performance. Our model incorporates 4 key components: spatial movement, pandemic spread, hospitalization dynamics, and warfare dynamics. By integrating these components, we offer a comprehensive representation of the interplay between pandemics and wars and their mutual impacts.
As an intervention strategy to reduce overall population deaths, we propose a DRL-based model that optimizes health care administration policy using in silico data generated by the proposed agent-based simulation framework. The inclusion of a dual health care system allows for the differentiation between civilian and military medical services, highlighting the unique challenges and administrative policies in both contexts.
As shown in Figure 5, the DRL model successfully learns an optimal administration policy, reducing the death toll by 70% compared to a naive policy in which patients are randomly allocated, provided hospitals are not at full capacity. Moreover, the evaluation of the DRL model compared to other baseline and single-usage health care policies, as shown in Figure 6, highlights its effectiveness. The DRL policy optimized for overall population health demonstrated statistically significant improvements over other strategies. Specifically, the allocation strategy focused solely on soldiers led to fewer deaths than the one focused solely on civilians, but the combined optimization yielded the best overall results. These findings underline the importance of accounting for both military and civilian health in designing administrative policies, as neglecting 1 population segment leads to suboptimal outcomes. In addition, focusing on only 1 population segment can raise sociological and ethical challenges in practice.
Further insights into the dynamics of death and injury under the optimal DRL administration policy are provided in Figure 7. The time-dependent behavior of wounded soldiers and overall deaths reveals non-linear and phase-dependent dynamics. For example, the sharp decline in wounded soldiers around the 185th simulation step coincides with the end of the war, whereas the sigmoid-like increase in total deaths reflects the dual impact of injuries from war and deaths from the pandemic. The delay in responding to pandemic needs due to the saturation of the health care system by injured soldiers highlights a critical trade-off between war-related injuries and pandemic management. To this end, Figure 8 illustrates the relationship between the optimization event horizon and the performance of the DRL model. The results indicate that shorter event optimization horizons, aligned with the average recovery rate, yield better outcomes. This suggests that the predictability of war dynamics, as described by the Lanchester model, allows the DRL model to focus on mitigating the more uncertain and impactful pandemic spread. Importantly, even with longer event optimization horizons, the model maintains a degree of robustness, though the influence of pandemic dynamics becomes more pronounced.
This study is not without limitations. First, the proposed model and simulation considered relatively short periods of time and therefore do not take into account population growth, reinfection, and multipathogen pandemics.Reference Alexi, Rosenfeld and Lazebnik 4 , Reference Khyar and Allali 41 , Reference Lopez and Rodo 57 , Reference Shami and Lazebnik 77 , Reference Turkyilmazoglu 84 However, because wars can be extended over many years,Reference Sullivan 81 this assumption may not hold. Future work could relax these assumptions to improve the model’s robustness. In addition, the model assumes simplified military contact patterns and neglects the complexity of modern warfare. Although the current approach assumes a historical perspective in which battles occur far from civilian zones,Reference Thorne 83 extending the model to include civilian casualties common in contemporary conflictsReference Khorram-Manesh, Burkle and Goniewicz 40 would improve its applicability. Furthermore, the model assumes isolation, ignoring international interdependencies during wars and pandemics.Reference Aguas, White and Hupert 2 Lastly, the simulation only considers 1 side of the war, neglecting the asymmetric nature of conflicts and assuming a static opposing force. Future studies should incorporate dynamic interactions between opposing forces to better capture real-world complexitiesReference Sabin 73 and consider more refined socio-demographic grouping for the administration policies.Reference Savchenko, Rosenfeld and Bunimovich-Mendrazitsky 75 , Reference Yaniv-Rosenfeld, Savchenko and Elalouf 95 , Reference Yaniv-Rosenfeld, Savchenko and Netzer 96
Taken together, our model provides valuable insights into designing effective health care and pandemic policies during wartime. For government agencies, understanding the predicted impact of conflict on disease spread can help prioritize resource allocation and intervention strategies. For military organizations, the model’s predictions can guide operational planning to minimize health impacts on both civilian populations in war zones and military personnel. Further investigation is required to refine and extend the proposed model for greater practical utility.
Data availability statement
The code and data that have been used in this study are publicly available in this study’s GitHub repository: https://github.com/shuchaa/pandemic_during_war_abs.
Acknowledgments
The authors wish to thank Noga Givon-Lavi for the clinical-epidemiological consulting.
Author contribution
Adi Shuchami: Conceptualization, Formal Analysis, Investigation, Software, Visualization, Writing – original draft.
Teddy Lazebnik: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Supervision, Visualization, Writing – original draft, Writing – review & editing.
Funding statement
This study received no funding.
Competing interests
The authors declare none.
