Distribution equality as an optimal epidemic mitigation strategy

Upon the development of a therapeutic, a successful response to a global pandemic relies on efficient worldwide distribution, a process constrained by our global shipping network. Most existing strategies seek to maximize the outflow of the therapeutics, hence optimizing for rapid dissemination. Here we find that this intuitive approach is, in fact, counterproductive. The reason is that by focusing strictly on the quantity of disseminated therapeutics, these strategies disregard the way in which this quantity distributes across destinations. Most crucially—they overlook the interplay of the therapeutic spreading patterns with those of the pathogens. This results in a discrepancy between supply and demand, that prohibits efficient mitigation even under optimal conditions of superfluous flow. To solve this, we design a dissemination strategy that naturally follows the predicted spreading patterns of the pathogens, optimizing not just for supply volume, but also for its congruency with the anticipated demand. Specifically, we show that epidemics spread relatively uniformly across all destinations, prompting us to introduce an equality constraint into our dissemination that prioritizes supply homogeneity. This strategy may, at times, slow down the supply rate in certain locations, however, thanks to its egalitarian nature, which mimics the flow of the pathogens, it provides a dramatic leap in overall mitigation efficiency, potentially saving more lives with orders of magnitude less resources.

1 Modeling the SIR epidemic spread and its mitigation 1

.1 SIR model
To model epidemic spreading via air travel 1 we used a local susceptible-infected-recovered (SIR) model with diffusive coupling, following the formulation presented in Ref. 2 . In this framework each node n (n = 1, . . . , N ) represents a local population of M n individuals, of which S n (t) are in the susceptible state S, I n (t) are infected (I) and R n (t) are recovered (R), hence S n (t) + I n (t) + R n (t) = M n for all t. Within each node, we assume a well-mixed population that locally follows SIR dynamics, [3][4][5][6] namely where α and β are the infection and recovery rates, respectively. The coupling between two populations n and m is mediated by the flux of incoming/outgoing travelers via air-traffic F nm , quantifying the number of individuals flying from m to n per day, namely F nm = F n←m .
This form of epidemic spreading gives rise to the following dynamic equations 2,7 dS n dt = −α S n (t)I n (t) M n + w nm R m (t) − w mn R n (t) . (1.4) The first term/s on the r.h.s. capture the processes of infection -proportional to the product of susceptible and infected individuals -and recovery, proportional to I n (t). The summation terms describe the diffusion of S, I or R individuals between local populations, where is the per-capita flux from m to n, hence, e.g., w nm S m (t) is the volume of susceptible passengers leaving m and entering n per day.
Finally, we introduce an invasion threshold ε, which activates the local SIR dynamics only if the infected population rises above an ε fraction of the local population. We apply this by adding a sigmoidal function 2 σ(x) = (x/ε) h 1 + (x/ε) h (1.6) to the relevant equation terms, providing dS n dt = −α S n (t)I n (t) M n σ I n (t) M n + N m=1 w nm S m (t) − w mn S n (t) (1.7) dI n dt = α S n (t)I n (t) M n σ I n (t) M n − βI n (t) + N m=1 w nm I m (t) − w mn I n (t) (1.8) dR n dt = βI n (t) + N m=1 w nm R m (t) − w mn R n (t) . (1.9) Hence the process of infection is locally initiated (i.e. turned on) only when I n (t)/M n exceeds ε, capturing an invasion of n.
Invasion threshold. The rational behind this threshold is to correct for an artifact that is otherwise unavoidable in our differential equation-based modeling. To understand this, let us first consider a stochastic realization of an SIR outbreak. It is seeded with a single infected individual, which then transmits the pathogen to its network neighbors. Quite often, due to the random nature of the spreading process, this initial seed may recover before it has the opportunity to transmit. This is the case even if the disease is in the pandemic phase with R 0 > 1, as long as the number of infection individuals -the seed -is small. Indeed, the sharp distinction between R 0 above or below unity is based on the mean-field analysis, relevant when the system is seeded with a macroscopic number of infections. Under such conditions the law of large numbers ensures a predictable decay (R 0 < 1) or proliferation (R 0 ≥ 1). With a small number of infections, however, the disease may taper of, simply by chance, before it has the opportunity to invade the network.
This realistic behavior, however, is deeply ingrained in the stochastic nature of real pandemics, especially at their early stages, but overlooked in the deterministic setting of our SIR equations. Indeed, under the differential equation-based modeling, even an infinitesimal seed δI n → 0 is guaranteed to instigate a local pandemic spread with 100% certainty. Once, however, the seed grows sufficiently, such that I n is no longer of order one or less, the mean-field solution kicks-in, and a stochastic elimination of the intial outbreak is no longer probable. This is precisely the role of σ(x), which only activates out deterministic spread when I n (t) M n > ε. (1.10) This rationale provides clear guidelines for selecting ε. Assuming a total of N destinations, that together account for the global population Ω, we have the average population being ⟨M ⟩ = Ω/N . Seeking I n (t) ≳ 1, and substituting M n in (1.10) by the average ⟨M ⟩, we arrive at ensuring an invasion of at least few individuals before the outbreak successfully penetrates a node.

SIR model with treatment
We now consider drug distribution, which divides the infected population into two separate groups, treated, who received the drug Q, and untreated, who have not yet acquired it, hence namely the infected population is constructed from the unity of treated and untreated infected individuals. We assume that uninfected individuals do not seek or receive treatment, i.e. no preventive measures, and hence the susceptible state S remains undivided. Similarly, once recovered the treated and untreated individuals reach the same state -which is non susceptible and non infectious -hence we also do not distinguish between treated and untreated individuals within the R population. Consequently, the possible transitions in each node are now driven by where ρ is the drug consumption rate, β is the recovery rate from the disease untreated, as appears in (1.1), and ζ > β is the recovery rate under treatment.
We arrive at the following four equations is the total infected population -treated and untreated -in n.
The consumption rate ρ depends on the drug availability in n, being zero in the limit of extreme shortage and approaching unity if there is sufficient availability to meet the demand. Hence we denote the total number of available drug doses in n by Q n (t), and its per-capita availability by , (1.19) in which the denominator represents the fraction of n's population that at time t continues to seek treatment. We then define the rate ρ as i.e. Q n (t) < I U n (t), while avoiding excess consumption when Q n (t) ≥ I U n (t), a case in which n has more available doses than required for its currently untreated population.
To close the set of equations (1.14) -(1.17) we must track the local availability of drugs Q n (t) at all locations. This is determined by our dissemination strategy, either Max-flow or Egalitarian, as detailed in Supplementary Section 2. Once selected, the dissemination strategy provides us with the daily influx L n← (t) of drugs into n at each day t (doses per unit time). Starting from t ≥ t R , our response time, we write where the first term represents the incoming drugs, as dictated by our dissemination scheme, and the second term captures their consumption. The Heavyside step-function ensures dissemination only begins at t = t R .

Normalized equations
Next we rewrite Eqs. (1.14) -(1.17) for the normalized populations s n (t) = S n (t)/M n , j U n (t) = I U n (t)/M n , j T n (t) = I T n (t)/M n and r n (t) = R n (t)/M n , obtaining where j n = j U n (t) + j T n (t). Assuming a negligible fraction of one-directional trips, i.e. that immigration accounts for a marginal part of the overall international mobility, we can write F nm = F mn , namely that, on average, the number of passengers flying daily from m to n is the same as those flying from n to m (see Supplementary Section 3). This enables us, using (1.5) to write is the normalized human flux matrix.
Finally, the normalized Eq. (1.21), takes the form where q n (t) = Q n (t)/M n is the normalized drug availability, set to q n (t) = 1 when there are sufficient doses for the entire population M n , and is the normalized influx, measuring incoming doses as a fraction of n's potential demand. Note that φ n in (1.19) and hence ρ in (1.20) remain unchanged following this normalization, i.e. φ n = Q n (t)/I U n (t) = q n (t)/j U n (t).

The human flux matrix
To construct A nm from data we denote the mobility rate by 2 in which the numerator quantifies the number of individuals flying per day and the denominator equals to the total population, therefore providing the daily fraction of individuals who seek air travel. This allows us to estimate a local population M n by 2 assuming that the total number of people departing from n, i.e. N k=1 F kn , represents a µ fraction of n's total population, M n . We can now write which allows us to extract A nm directly from the mobility data (F nm , µ).

Summary -SIR dynamics
Normalized equations for epidemic spreading with treatment: The relevant parameters: Passenger fluxes are characterized by F nm and µ, evaluated from the flight data • Disease/therapeutic parameters α, β and ζ are selected to represent different epidemiological scenarios, from benign to severe.
2 Dissemination of therapeutics

Overview
To examine our dissemination strategies we used the commodity flow framework 8 to obtain the optimal shipping strategy. We consider the air-traffic network as our underlying flow network, upon which the single source node s must ship the therapeutic in sufficient quantities, to meet all initial demands D n (0) of all target nodes n = 1, . . . , N . Each link is characterized by a non-negative carrying capacity B nm , capturing the maximal volume of the therapeutic that can be transferred from m to n through that link per unit time (day). Our goal is to meet the demand of all nodes in an optimal fashion, i.e. optimizing for Max-flow or Egalitarian, within the shipping constraints imposed by B nm .

Optimization
At t = 0 we select a source node s, and set the initial demands of all nodes to D n (0), n = 1, . . . , N , capturing the number of individual therapeutic (drug or vaccine) doses required at n. Different strategies for setting D n (0) are discussed bellow. Denoting the volume of doses that are actually shipped from m to n on day t by L nm (t), we seek to maximize the total number of doses shipped each day from s. Under the egalitarian strategy, we also wish for these doses to be disseminated as evenly as possible across all nodes. Together, this translates to is the net outgoing flux of drug doses from n at time t, and is the incoming flux to n, as appears in Eq. (1.21).
The second term in (2.1) captures the level of homogeneity in each dissemination step t. First, we denote by λ n (t) = L n← (t)/D n (t) the fraction of n's demand that was supplied to n at time t. We then use this to obtain which is upper-bounded by the least supplied node at time-step t. Under a perfectly ho-mogeneous spread of the therapeutic we have a uniform λ n (t), i.e. λ n (t) = λ(t) for all n, and hence L s→ (t) = N n=1 λ(t)D n (t). Under these conditions the second term of (2.1) has a comparable contribution to the optimization as that of the first. In contrast, an uneven dissemination will lead to some λ n (t) being smaller, bounding the value of λ(t), and consequently decreasing the target function (2.1). For example, a sequential supply pattern, in which nodes are supplied serially, is characterized by one or more destinations with λ n (t) = 0 at each shipping instance. Indeed, such dissemination scheme first rapidly fills up the demand of few nodes, then turns to fill up that of the rest. Such supply strategy, while potentially fast, is highly uneven and has λ(t) = 0 at most time-steps. Equation (2.1) penalizes such inequality by having the second term vanish.
Therefore, the first term in (2.1) optimizes for Max-flow, while the second term enhances Egalitarian dissemination. The relative role of these two competing target functions is governed by the homogeneity coefficient H, whose value determines whether we prioritize flow volume (H → 0) or flow equality (H ≈ 1).
Following each day, the demands are updated as subtracting the incoming supply from the previous day's demand, thus obtaining the remaining demand at t + 1. The maximization of (2.1) is subject to three constraints that must be satisfied for all t: The first constraint (2.6) ensures that all daily shipments from m to n do not exceed the carrying capacity B nm ; the second constraint (2.7), with D n (t) taken from (2.5), restricts nodes from accumulating therapeutics in excess of their remaining demand, therefore uperbounding the total incoming doses at day t, by the current demand D n (t). Constraint (2.8) is only relevant in case H ̸ = 0 in (2.1). Its goal, in accordance with Eq. (2.4), is to ensure all nodes receive at the least a portion λ(t) of their demand at every time-step, avoiding a supply pattern in which few nodes receive their demand early on, while the rest are only supplied at a later stage. Under H = 0, the value of λ(t) is irrelevant in the optimization, and therefore the optimal solution may result in λ(t) = 0, allowing a selective, and hence unequal, supply pattern.
Carrying capacity. To construct the carrying capacities B nm we used the empirical passenger fluxes F nm , as obtained from the air travel data 9 . We assume that the daily volume of cargo shipped across the m → n route is proportional to the number of passengers F nm for that route 10 . This allows us to write where κ is the cargo conversion coefficient, capturing the number of doses that can be shipped as cargo per each traveling passenger. Our results are not sensitive to the specific value of κ, as any re-scaling of B nm (capturing doses/day) is equivalent to an arbitrary re-scaling of the time units, having no bearing on our essential findings. Moreover, it is also possible that shipping takes place on different routes altogether, and hence the shipping network may be distinct from the travel network. Fortunately, our results do not depend on the specific network structure, as we show, using an arbitrary random network in Supplementary Section 5.

Normalized flow equations
which we further express as are the per-capita (in n) incoming rates from m to n. In (2.12) and (2.13) the population ratio M m /M n can be evaluated from the data based on Eq. (1.36), providing expressed in terms of the empirically available human flux matrix F nm . These normalized rates, l n← (t) and l n→ (t), are subject to the shipping constraints, which can be rewritten by substituting (2.14) and (1.34) into (2.6) -(2.8), providing In these normalized expressions the carrying capacities are and the per-capita demands are Finally, dividing both sides of Eq. (2.5) by M n , we obtain and the optimization of Eq. (2.1) becomes Max-flow. To optimize for max-flow we set H = 0 in (2.22). Under these conditions the optimization scheme seeks to only maximize the net outflux from s, l s→ (t), namely obtain a shipping sequence in which the volume of therapeutics introduced at each time-step into the network is as high as possible. Such sequence may be, and in fact, typically is, highly uneven. For example, supplying the entire demand of few nodes, then proceeding to supply that of the remaining nodes. Such supply pattern, while globally efficient, may be very inefficient for the supply-delayed nodes.
Egalitarian. Setting H ≈ 1 we now favor supply sequences that also satisfy constraint (2.18). Such optimization penalizes unequal supply patterns, since the value of λ(t) is determined, at each iteration, by the minimally supplied node in the network, as expressed in (2.4). Hence the optimal dissemination seeks not just a high volume of drugs shipped daily, but also an even spread of that volume, for which λ(t) is maximal.
Normalized carrying capacities. To construct b nm from data we take B nm from (2.9) and M n from (1.36), which in (2.19) provides us with (2.23) Hence, b nm depends on the passenger data, as captured by the fluxes F nm and by the mobility rate µ -parameters that have been already introduced in the construction of the SIR dynamics. It requires, in addition, the conversion ratio κ, which controls the absolute numbers of doses that can be carried by each shipping route.
Normalized capacity C s . The overall shipping capacity of the air transportation network is governed by the conversion ratio κ, capturing the volume of cargo (in doses) that can be shipped across each route per human passenger. In our presentation, however, we used the normalized C s ∝ κ, to describe the network's capacity in more natural units. Consider the maximal number of doses that can potentially be shipped from the source in each day. Following the constraint (2.6) we can express this as capturing the fact that even if our optimization achieves ideal shipping rates, the outgoing flux of drugs from s is capped by the carrying capacity of all of s's outbound flight routes. Dividing by the global population, Ω = N n=1 M n , we obtain the daily per-capita shipping potential of s, namely C s = max(L s→ (t))/Ω. We rewrite this as where we used (2.9) and (1.36) to express the numerator and the denominator, respectively.
The normalized capacity (2.25) is, indeed, proportional to κ, with a multiplicative factor that is fully determined by the available mobility data (F nm , µ). It is gauged to follow C s = 1 in case the network can, under ideal conditions, achieve an outflux of order Ω, i.e. disseminate sufficient doses to supply the global demand per day. Note, that this interpretation of C s is independent of s, hence a capacity of, e.g., C s = 1 has the same meaning regardless of the identity of the source node. For example, consider two source nodes s and s ′ whose absolute shipping capabilities are N m=1 B ms = 2 N m=1 B ms ′ , namely s has twice the shipping capabilities of s ′ . Under any specific choice of κ we have C s = 2C s ′ , capturing the fact that, indeed, under similar conditions s will achieve twice the dissemination efficiency as compared to s ′ . Conversely, if we wish to compare performance under a preset C s , as we do in the main text, we must set κ ′ = 2κ, compensating for the difference between the two distinct source nodes. Hence C s allows us to compare different mitigation scenarios without being sensitive to the specific characteristics of each realization, for instance the arbitrary selection of s.

Setting the local demands
We considered three different strategies for assigning the initial demands d n (0): Population based. We set the demand in each location based on that location's population, i.e. D n (0) = M n . Under the normalized flow dynamics this maps to d n (0) = 1 for all nodes. In our formulation, since M n in (1.36) is proportional each node's weighted degree, this is equivalent to a degree-based supply pattern. Demand based. In the demand based strategy we used the projected infection rate at each destination to evaluate its anticipated demand. Hence, we first simulated the unmitigated epidemic via Eqs. (1.40) -(1.43) with q n (t) = 0. Then, extracting r n (t → ∞), we obtain the fraction of each population that will, at some point, become infected and require treatment.
To meet this demand level we set d n (0) = r n (t → ∞), hence we supply each node based on its predicted demand.
Urgency based. Here, we implemented a dynamic strategy, in which the demands are updated in real-time, based on the states of infection (j n (t)) and drug availability (q n (t)) at each node. Therefore, we replace Eq. (2.5) with a dynamic update rule, following d n (t + 1) = j n (t) − q n (t) − l n← (t). (2.26) Here the demand is updated at each time-step (day) based on the supply gaps observed in the previous time-step. This is captured by the difference between j n (t), capturing the demand level at time t, and the actual supply, which equals to the previous day supply l n← (t) added to the already available drugs q n (t). The value of j n (t) and q n (t) are extracted from the SIR equations at the beginning of each day. In Eq. (2.26) time is measured starting from t R , hence the initial demands are determined by the state of the epidemic at t R .
Results obtained from Max-flow and Egalitarian with all three strategies are presented in Fig. S1. In all cases egalitarian (green) outperforms max-flow (yellow), as expected. Under the urgency-based scheme (Fig. S1c) the margin between max-flow and egalitarian becomes small. This is because this version of max-flow, by design, coordinates supply with demand. Indeed, each day's supply is set according to the demand left from the previous day at the relevant destination. Therefore, we arrive at a supply pattern that naturally mimics the spread of the pathogens, in a sense, rendering max-flow to deliver therapeutics in a similar fashion to that of egalitarian. In other words, the urgency-based scheme naturally leads to egalitarian supply, as it emulates the homogeneous spreading patterns of the pathogens. Yet, while the urgency-based strategy requires real-time update and tracking, egalitarian constructs the required homogeneous supply patterns quite simply from the start.

Summary -Flow dynamics
Normalized equations and constraints for flow optimization: Maximize l s→ (t) + Hλ(t) . where the preset demands 0 ≤ d n (0) ≤ 1 are updated via To obtain the optimal solution for (2.27) we use an iterative approach (Fig. 1, main text): at each time step t we use linear optimization to maximize l s→ (t) under the three constraints. This provides us with day t's optimized shipping strategy, detailing the daily fluxes l nm (t) through all routes m → n. With these fluxes we calculate the net influx into each node l n← (t), and update the demands d n (t + 1). We then use these updated demands to obtain the shipping strategy for the next day t + 1, repeating the process until all demands vanish. The supply time T Supp,n of node n is reached when d n (T Supp,n ) = 0 for the first time.

Mobility networks and epidemic parameters
Passenger mobility network. We used the Global Mobility Network dataset 9 to extract the empirical fluxes F nm . The dataset lists the daily travel of ∼ 8.91 × 10 6 passengers per day over the course of 365 days, i.e. i = 1, . . . , 365, between N = 1, 292 airports, linked through 38, 377 directional air-routes. The mean flux from m to n is thus captured by where F {i} n←m is the number of passengers flying from m to n in day i. In our analysis we assume an almost symmetric flux matrix, as expressed in Eq. (1.27). Indeed we find that practically all links are bi-directional, and that in the vast majority of existing links we have F n←m ≈ F n→m , a similar passenger flux in both directions. Therefore, to construct F nm we used the symmetrized 2 Disease parameters. The disease parameters in (1.40) -(1.44) were set to β = 0.285 day −1 and α = 1.5β, 2β and 10β, providing a reproduction rate of R 0 = α/β = 1.5, 2, 10, as appears in the main text, covering different levels of contagion 3 . The treated population recovery rate is set to ζ = 1 day −1 . For the penetration threshold σ(j n ) in (1.6) we used Eq. (1.11) with N ∼ 10 3 and Ω ∼ 10 9 to set ε = 10 −6 ; the saturation coefficient was set to h = 8. The response time t R was set to be at the onset of the pandemic spread, evaluated as the time t when the infection levels have reached 10% of the projected peak infection (Fig. S2). The above represents our default parameters, used, e.g., to produce Fig. 2c,d of the main text; in other figures we consistently vary some of these parameters, such as t R or R 0 as reported in the main text. Excluding the varied parameter in each panel, e.g., t R in Fig. 2i, the remaining parameters continue to be fixed at the above default values.

Numerical analysis
We used a 4th order Runge-Katta stepper to iteratively solve the SIR dynamics of Eqs.

spread. To numerically realize this limit we implemented the termination condition
where t i is the time stamp of the ith Runge-Kutta step and ∆t i = t i − t i−1 . This condition tracks the ratio between the numerical derivativeṡ n = s n (t i )−s n (t i−1 )/∆t i and s n (t i ), guaranteeing that its maximum value over all nodes is smaller than the pre-defined termination parameter ξ. In our simulations we set ξ = 10 −3 .
To model our mitigation we used linear optimization to maximize (2.27) in Eqs. (2.27) -(2.32). The optimization yields a supply sequence in the form of l nm (t), capturing the normalized daily therapeutic fluxes between all pairs of nodes m → n. We then use (2.32) to extract the daily influx l n← (t), quantifying the number of doses per capita that enter node n at each day t. The obtained l n← (t) were then used as input in Eq. (1.44) to track the local drug availability in all nodes. While the Runge-Kutta stepper, implemented to solve the epidemic dynamics, incorporates variable time steps, typically much shorter than a single day, the optimization scheme provides l n← (t) at a daily resolution. Therefore, to couple the commodity flow equations with the SIR dynamics we first applied the flow dynamics to extract l n← (t). We then advanced the SIR equations day by day, updating q n (t) at the beginning of each day as q n (t + 1) = q n (t) + l n← (t), (3.4) after which q n (t) continues to evolve via (1.44). To be specific we used the following procedure: Set the response time t R (days). As no drugs were supplied prior to this time point we set l n← (t) = 0 for t = 0, . . . , t R .
Use the commodity flow optimization to obtain l n← (t) for t = t R + 1, . . . , T , where T is the overall supply time. The result is the inflow of drugs into n at each day t.
Use variable time-step Runge-Katta to solve the SIR equations day by day. For example, at t = 0 we advance the equations until reaching the end of the first day at t = 1. At t = 1 we update q n (t) via (3.4), and set the resulting q n (1) as the new initial condition for the next day. We then continue to solve the SIR dynamics with initial condition q n (t) from t = 1 to t = 2. We continue the process until the epidemic reaches saturation, i.e. satisfies (3.3).
This process captures the continuous nature of the epidemic spread (variable time-step Runge-Kutta), and its coupling with the discrete daily bursts of the supply intervals. Note that since l n← (t) = 0 for t < t R the day by day advancement of step 3 above only comes into effect when t > t R .
This coupling of the continuous-time SIR equations with the discrete-time commodity-flow dynamics becomes compatible in the limit where the discrete time steps are small compared to the SIR natural timescales. Under this limit the state of the pandemic does not change significantly between supply instances, and the discreteness is smoothed out. Here, as the pandemic timescales are of order of ∼ 10 2 days, our day-by-day tracking of the supply, indeed, satisfies this condition.

Heterogeneity in network-based distribution
Consider the flow of a unit therapeutic from a single source node s, as it traverses all network pathways. With each step along its trajectory, it becomes diluted, first splitting among s's nearest neighbors, then its second neighbors and so on. In a random network, due to the locally tree-like structure, the number of nodes at distance r from the source grows exponentially as 11 where λ depends on the specific network topology. Hence, the therapeutic flux l s→ (t), exiting s daily, is dispersed among an exponentially growing number of nodes as it penetrates the network. For a typical node at distance r from s we, therefore, expect a supply rate of 12 S(r) = l s→ (t)e −λr (4.2) an exponential decay that ensures a conservation of all doses flowing from s.
We can now use (4.2) to obtain the density function P (S) to observe S ∈ (S, S + dS) by transforming from the random variable S to the variable r. We write Taking r from (4.2) we have r = −λ −1 ln S and hence | dr/ dS| ∼ S −1 . Using (4.1) and (4.2) we can also write P (r) ∼ S −1 , hence, together we arrive at where ν = 2, precisely the prediction of Fig. 3d of the main text.
This result is obtained under the configuration model ensemble, i.e. random network/weights. More generally, e.g., in the presence of degree/weight correlations, the exponent ν can take different values. Similarly, if drugs are consumed as they flow through the network, a process neglected in the above derivation, the decay in (4.2) may become more rapid, again, impacting the exponent ν. Still, the power-law form of (4.4) will remain unaffected, albeit having ν ̸ = 2. Indeed, as long as the network exhibits an exponential growth of the form (4.1), and the downstream dilution exhibits an exponential decay of the form (4.2), we expect a scale-free P (S). To observe this, consider, for example, that the two exponential functions exhibit different growth/decay rates, i.e. P (r) ∼ e λ 1 r (4.5) Under these conditions we continue to have | dr/ dS| ∼ S −1 , however P (r) now follows P (r) ∼ S −λ 1 /λ 2 . Using these expressions in (4.3), we find that the scaling P (S) ∼ S −ν remains, only that now the scaling is which may assume values different than 2. In case λ 1 = λ 2 , Eq. (4.7) converges to our original derivation and provides ν = 2. In Fig. 5d of the main text, we found that the empirical supply rates of the COVID-19 vaccine dissemination followed (4.4) with ν = 1.7, consistent with (4.7) under λ 1 /λ 2 = 0.7.

Additional tests
Our analysis of the two distributions P (I) and P (S) is rooted in the interplay between the network structure, i.e. the exponential growth of (4.1) and the spreading dynamics. For the pathogen, this exponential growth leads to the extremely short pathways, that help the disease spread to the majority of nodes concurrently, resulting in the bounded P (I). For the therapeutics it leads to the exponential dilution of (4.2), which, in turn, renders P (S) scale-free. The crucial point is, that (4.1) is a quite generic characteristic of, practically all random networks, and therefore, as we explain below in detail, we expect our analysis to be largely independent of the specific network structure, or of the arbitrary selection of the source node.
To examine this more deeply, we consider the impact of the source node selection. In a heterogeneous network, e.g., scale-free, different nodes may, in principle, have a distinct local environment, and hence one might expect that they would lead to a different dissemination outcome. Still, while such differences are local, capturing a node's immediate surrounding, they do not lead to different global dynamics. This is because the exponential growth rate λ in (4.1) captures a global network characteristic, which is largely independent of the source node. To understand this consider two highly distinct nodes: node 1 a hub, and node 2, a peripheral node. At r = 1, these two nodes interact with a profoundly different environment, with 1 having many neighbors and 2 having just a few. However, as we advance further at r = 2, 3, . . . the more distant neighbors of both 1 and 2 become statistically indistinguishable. 11 Therefore, the exponential growth rapidly converges to (4.1), with the network's global growth rate λ. Consequently, the global spreading patterns of the pathogens (P (I)) and the therapeutics (P (S)) are independent of (i) the specific network structure; (ii) the source of the disease outbreak; and (iii) the source of the therapeutic distribution. Below we examine each of these predictions (i) -(iii).
Network structure. To observe the fact that network structure has little impact on the supply/demand patterns we constructed an Erdős-Rényi network of N = 1, 000 nodes with average degree ⟨k⟩ = 15, and link weights extracted from a normal distribution N (µ, σ 2 ), setting µ = 1, σ 2 = 0.1 (and discarding negative weights). This represents a fundamentally different network from the air-transportation network used in our main analysis, which has a scale-free topology and extremely heterogeneous weights. Still, as Fig. S3 clearly shows the pathogen/therapeutic distributions continue to follow our predicted patterns: the pathogens exhibit a homogeneous spread, and hence a concurrent global demand, while the therapeutics follow the predicted power-law (ν = 2), depicting a highly unbalanced supply.
Outbreak/dissemination source. To complete our examination of the max-flow vs. egalitarian dissemination, we repeated our analysis this time selecting a set of 6 random assignments for the outbreak node o and the therapeutic source node s (Fig. S4). In all cases we examine the mitigation efficiency E vs. capacity C, a là Fig. 2h of the main text. As predicted, we observe no discernible difference between the resulting plots, as -indeed -our analysis is independent of the selection of source/outbreak nodes.
Multiple sources. In a real-world setting, the therapeutic may be distributed from more than one source s. Using multiple sources, can, quite naturally, enhance the performance of both strategies -max-flow and egalitarian. However, as long as the number of sources n is small, i.e. n ≪ N , the egalitarian advantage is guaranteed. This is because, in a large network, each source continues to be surrounded by a tree-like expansion of the form (4.1), leading to the predicted supply/demand discrepancy. Of course, when the number of sources is increased, as we approach n ≲ N , we have each node being in close proximity to one of the sources. Under these conditions, all nodes benefit from a similar supply rate -indeed they all have a supplier in their neighborhood -and hence they are all equally treated. As a consequence, for all purposes, such conditions converge to our proposed egalitarian dissemination, and the max-flow/egalitarian distinction is rendered moot. In Fig. S5 we examine the impact of 2 or 3 distribution sources, finding that n has, indeed, little impact on the egalitarian advantage.

Vaccine-based mitigation
Egalitarian mitigation is designed around drug distribution, but continues to provide improved performance, albeit less dramatic, also in the case of vaccines. To examine this we implemented a susceptible-infected-vaccinated-recovered (SIVR) model, in which the disease is battled through vaccines rather than drugs. The resulting transitions are now and in which φ n (t) is the vaccine availability as a fraction of the potential demand, here captured by the susceptible population.
where s n , j n , v n and r n are the normalized S, I, V, R populations, and q n is the normalized vaccine availability. Simulating (5.6) -(5.10), in Fig. S6 we extracted the efficiency E vs. the shipping capacity C under both max-flow (yellow) and egalitarian (green). We find that egalitarian offers improved mitigation, however, the margin is smaller than that observed under drug distribution.

Simulating COVID-19
To track the COVID-19 scenario we added an additional compartments to our SIR modeling framework: the exposed individuals (E), who have contracted the virus, but are still not infectious. The resulting transitions are now where s n , e n , j n and r n are the normalized S, E, I, R populations, and q n is the normalized drug availability. The air-transportation fluxes A nm are extracted from our mobility data, via Eq. (1.39), namely allowing us to tune the volume of international travel by changing µ, in order to account for the observed reduction in air travel in response to the pandemic. In (6.4) and (6.5) we omitted the travel terms completely for infected individuals.
Evaluating the parameters. To construct a realistic simulation we collected data on the time-scales of the different transitions characterizing the COVID-19 disease cycle. 13 This provides us with γ = 1/3 days −1 and β = 0.1 days −1 , capturing the average 3 day Fitting the data to a shifted exponential of the form (6.9) we extract the growth rate χ n and the delay τ P n . (d) The probability density P (χ) for a random country to have growth rate χ n ∈ (χ, χ + dχ). (e) The probability density P (τ P ) for a random country to have delay τ P n ∈ (τ P , τ P + dτ P ). (f) The probability density P (I) for a random country to have infection rate I n ∈ (I, I + dI).
incubation period, and 10 day recovery time-line. To capture the reduction in air-travel we set µ = 0.85 × 10 −3 day −1 , a reduction of 50% as compared to normal circumstances.
The infection rate α cannot be directly measured, and we therefore evaluated it from the observed spreading dynamics. Hence, we collected data from the early stages of the pandemic in 41 different countries. 14 In each country we track the initial exponential rise in j n (t), prior to the instigation of lock-downs or quarantines (Fig. S7a-c, Supplementary Dataset 1). We then fit the observed data to a shifted exponential of the form j n (t) = j n (0)e χn(t−τ P n ) , (6.9) in which χ n is the observed growth rate of the infected population at n and τ P n is the pathogen delay, capturing the propagation time for the virus to reach n. The pre-factor j n (0) is extracted directly from the first data-point in each country's epidemic data. The distribution of χ n across all 41 monitored countries is shown in Fig. S7d, a bounded distribution with mean growth rate of χ = 0.25 days −1 , consistent with other similar analyses. 15 Next we implemented our SEIR disease cycle with the relevant parameters γ = 1/3, β = 0.1, under q n (t) = 0, i.e. no vaccine, to seek the infection rate α than best matches the observed χ = 0.25. This allows us to infer α = 0.62 days −1 , completing our empirical extraction of all parameters required for Eqs. (6.2) -(6.7).
Extracting P (I) and P (S). The parameter τ P n in Eq. (6.9) captures the propagation time of the pathogen from the outbreak source o to the destination n. Hence, by fitting (6.9) to the observed spreading dynamics, we can extract the time needed for the pathogen to reach each of our 41 tracked destinations. The rate by which the pathogens advance towards n is, therefore (6.10) providing an empirical evaluation of Eq. (8) of the main text. In Fig. S7e we show P (I) as obtained from our empirical pool of countries, finding that it is indeed bounded, indicating the homogeneous SARS-CoV-2 spreading patterns.
To obtain the vaccine supply patterns, we collected data 16 on the vaccination rate of 229 countries (Supplementary Dataset 2). We then measured the elapsed time τ S n from the initial vaccine dissemination (December 1,2020), until country n began its vaccination campaign. This helps us extract the supply rate of Eq. (9) in the main text via S n = 1 τ S n , (6.11) whose distribution is displayed in Fig. 5d of the main text.