Digitizable therapeutics for decentralized mitigation of global pandemics.

When confronted with a globally spreading epidemic, we seek efficient strategies for drug dissemination, creating a competition between supply and demand at a global scale. Propagating along similar networks, e.g., air-transportation, the spreading dynamics of the supply vs. the demand are, however, fundamentally different, with the pathogens driven by contagion dynamics, and the drugs by commodity flow. We show that these different dynamics lead to intrinsically distinct spreading patterns: while viruses spread homogeneously across all destinations, creating a concurrent global demand, commodity flow unavoidably leads to a highly uneven spread, in which selected nodes are rapidly supplied, while the majority remains deprived. Consequently, even under ideal conditions of extreme production and shipping capacities, due to the inherent heterogeneity of network-based commodity flow, efficient mitigation becomes practically unattainable, as homogeneous demand is met by highly heterogeneous supply. Therefore, we propose here a decentralized mitigation strategy, based on local production and dissemination of therapeutics, that, in effect, bypasses the existing distribution networks. Such decentralization is enabled thanks to the recent development of digitizable therapeutics, based on, e.g., short DNA sequences or printable chemical compounds, that can be distributed as digital sequence files and synthesized on location via DNA/3D printing technology. We test our decentralized mitigation under extremely challenging conditions, such as suppressed local production rates or low therapeutic efficacy, and find that thanks to its homogeneous nature, it consistently outperforms the centralized alternative, saving many more lives with significantly less resources.


.1 SIR model
To model epidemic spreading via air travel [1] we used a local susceptible-infectedrecovered (SIR) model with diffusive coupling, following 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. This form of epidemic spreading gives rise to the following dynamic equations [2,7] dS n dt = −α S n (t)I n (t) dR n dt = βI n (t) + N m=1 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 infection is locally initiated only when I n (t)/M n exceeds ε, an invasion of n.

SIR model with treatment
We now consider drug distribution, which divides the population into two separate groups, the treated population, who received the drug Q, and the untreated population who have not yet acquired it. This leads to six distinct states, characterizing each node's population: where ρ is the drug consumption rate, β is the recovery rate from the untreated disease, as appears in (1.1), and ζ > β is the recovery rate under treatment. The drug efficacy is 0 ≤ γ ≤ 1, capturing the probability of recovery (at rate ζ) following drug consumption, hence for each consumed dose, at a rate ρ, a fraction γ of the consumers transitions to the treated state, while a 1 − γ fraction remains untreated. In other words, the first four reactions represent the process of treatment, in which a single dose Q is consumed and an untreated individual transitions with probability γ to the treated state. For a fatal disease, the untreated population perishes, hence R U n (t), the number of individuals at the R U state in n represents n's deceased population, while R T n (t), counting individuals in the R T state, represents the number of saved individuals, who would, in the absence of the therapeutic, also perish. In (1.12) we assume preemptive drug dissemination, in which once the therapeutics become available, they are provided to the entire populationsusceptible or infected -such that a susceptible individual can undergo treatment directly upon infection, without the additional lag for acquiring the therapeutic.
We arrive at the following six equations where is the total infected population -treated and untreated -in n. Note that we have now eliminated the diffusion of the R U population in (1.17), as this state represents the individuals who have died from the disease, and are thus redacted from the traveling population.
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.20) in which the denominator represents the time-dependent demand in n. We then define the rate ρ as ρ(ϕ n ) =    ϕ n for ϕ n < 1 1 for ϕ n ≥ 1 = min{ϕ n , 1}, (1.21) which in Eqs. (1.13) -(1.16) ensures that drug doses are consumed as long as there is local demand, while avoiding excess consumption when Q n (t) > S U n (t) + I U n (t), in which case n has more available doses than required for its current untreated population.
To close the set of equations (1.13) -(1.18) we track the drug availability Q n (t) in each location n. Under centralized mitigation drug doses are manufactured and shipped from a single source node s at a rate C s (doses/day), providing where δ ns = 1 if n = s and δ ns = 0 otherwise. The delay function θ is the Heavyside step-function initiating drug production only at t ≥ t R , incorporating the response time required before beginning drug distribution. The diffusion rates B nm represent the fraction of drug-doses in m that ship daily to n; see Sec. 2.2. Under decentralized mitigation drugs are synthesized locally, with each node characterized by its individual production rate C n , providing Hence, Eqs. (1.13) -(1.18) with ρ taken from (1.21) capture the spread of the epidemic and its competition with the distributed therapy. The distinction between centralized and decentralized mitigation is provided by the Q n (t) equations, taken from (1.22) or (1.24), respectively.

Normalized equations
Next we rewrite Eqs. (1.13) -(1.18) for the normalized populations s U n (t) = S U n (t)/M n , s T n (t) = S T n (t)/M n , j U n (t) = I U n (t)/M n , j T n (t) = I T n (t)/M n , r U n (t) = R U n (t)/M n and r T n (t) = R T 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 , stating 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 2). This enables us, using (1.5) to write is the normalized human flux matrix.
Next we normalize Eqs. (1.22) and (1.24) by tracking the per-capita drug availability q n (t) = Q n (t)/M n , which approaches unity when in n there is sufficient dosage for the entire local population. For centralized mitigation we obtain (1.22) where ρ is taken from (1.21) and in (1.20) remains unchanged under the normalization. The matrix represents the incoming flux of therapeutics from m to n. We define the normalized central production rate as which approaches C s = 1 when s has the capacity to produce and ship enough doses for the entire global population each day. Under common conditions we expect C s < 1.
Using this normalization Eq. (1.39) becomes where is the global population. Note that we used the δ-function in the second term to substitute M n in the denominator of (1.39) by M s in (1.43). Our normalization introduces the prefactor Ω/M s , which quantifies the global population in units of s's local population. This factor arises because under centralized mitigation, to meet the global demand, the source s must export doses in quantities that are of order Ω/M s greater than its local population. Therefore, a normalized capacity of C s = 1 translates to a per-capita in s production of Ω/M s doses per day. Such level of production, which does not scale with M s , but rather with Ω is, in the majority of cases, significantly higher than practically attainable, hence we expect, in real scenarios, that C s < 1. Under decentralized mitigation, Eq. (1.24) provides where c n = C n /M n is the per-capita local synthesis rate. Here c n = 1 translates to a production and local distribution capacity of M n doses per day, namely n can provide therapy to its entire local population. As opposed to C s , this normalization scales with the local, rather than global, population, hence C s = 1 in the single source node s and c n = 1 over all n = 1, . . . , N , both capture the case where the global demand can be met within one day -in the former in a centralized fashion and in the latter, decentralized across all N nodes.

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 assuming that all people departing from n, i.e. N k=1 F kn , represent a µ fraction of n's total population, M n . We can now write which allows us to extract A nm directly from mobility data (F nm , µ). A similar derivation, using Eq. (1.47), helps us express the ratio Ω/M s in (1.43) as The parameter κ s captures s's fraction of the global population. For a typical source node s, it scales as the number of nodes in the air-traffic network.

Human mobility network and epidemic parameters
We used the Global Mobility Network dataset [8] 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.31). 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 a network comprising 19, 614 bi-directional links, that represents a symmetric approximation of the empirical fluxes. Evaluating the local populations M n as in (1.47) we obtained A nm via (1.50). Following Ref. [2], we used µ = 1.7 × 10 −3 day −1 for the mobility rate. The disease parameters in (1.53) -(1.58) were set to α = 2 day −1 and β = 0.2 day −1 , providing a reproduction rate of R 0 = α/β = 10, and hence ensuring a highly contagious pandemic, which, untreated, will be guaranteed to spread globally [3]. Other, less extreme scenarios, are examined in Sec. 4. The treated population recovery rate is set to ζ = 1 day −1 , and the efficacy is set to γ = 1. For the penetration threshold in σ(j n ) we used ε = 10 −6 and h = 8 in (1.6). The term Ω/M s = κ s in (1.59) was evaluated using (1.51). Finally, 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. S1). For the parameters listed above, this provides t R = 12 days. As our source node we selected s=ITM (Itami International Airport, Osaka, Japan), testing the effect under other sources in Section 4.
The capacities C s and c n were set variably to test the performance of each method, centralized vs. decentralized with varying global/local capacities. In the case of decentralized mitigation the capacities are distributed, potentially taking different values in different locations n. Therefore, if e.g., we set the mean capacities to C s = c n = C, then C s = C, but c n ∼ N (C, σ 2 ) is a random variable extracted from a normal distribution with mean C and standard deviation σ = 0.1C.
The above represents our default parameters, used, e.g., to produce Fig. 1 of the main text; in Figs. 2 -3 we consistently vary some of these parameters, such as t R or γ as reported in the main text. Excluding the varied parameter in each panel, e.g., t R in Fig. 2a, the remaining parameters continue to be fixed at the above default values. The setting used to produce Fig. 4 of the main text are explicated in Sec. 3. In Sec. 4 we examine the behavior of our system under different disease parameter sets, varying the reproduction rate R 0 = α/β.  Figure S1: Selecting the response time. To select t R we tracked the global level of infection, I(t) = N n=1 I n (t), vs. t, and set t R to be the time where I(t) has reached an η fraction of its peak value, namely I(t R ) = ηI(T Peak ). Here we set η = 0.1.
In our simulations we track the epidemic spread using a 4th order Runge-Kutta stepper to iteratively solve Eqs. (1.53) -(1.60), with the relevant parameters as described above.

Drug distribution networks
The networks Z nm and B mn in (1.59) represent the fluxes of therapeutic doses to and from n respectively. These fluxes depend on the distribution strategy, which in turn impacts the efficiency of the drug dissemination. We consider three different constructions of these networks -two that build on the existing air-routes, and a third that achieves the highest efficiency by exerting extensive control over the international aviation network.
Diffusive distribution (Network 1). The results presented in the main text were obtained by setting This represents diffusive distribution, in which the drug spreads along similar routes as the human travelers. The rationale is that in each node m the q m (t) available doses are split among each of m's outgoing destinations. The fraction of these doses to reach n is equal to the fraction of passengers departing from m whose destination is n, namely we calculate the flux from m to n (numerator) as a fraction of the total flux outgoing from m (denominator). This natural distribution scheme requires no intervention in the aviation networks, as the busy routes, mobilizing many human travelers, also naturally ship a larger fraction of commodities, in this case -drug doses. It also leads, on average, to a fair distribution, as airports with large incoming fluxes, which are typically highly exposed to the spreading pathogens, also receive increased fractions of the therapy. Under this construction, the networks Z nm = B mn are symmetric, as they can both be expressed, using Eq. (1.50), as capturing an effective travel network in which 100% of passengers, in this case the available drug doses, travel each day (equivalent to setting µ = 1 in (1.50), i.e. a mobility rate of 100%).
Directed distribution (Network 2). As an alternative to the above network we consider a direct flight network, in which the therapeutic is shipped in directed routes from s to all destinations. First we consider the binary air-route network G nm , which satisfies capturing all existing air-routes, regardless of their volume of traffic. As the therapeutic is shipped from s to all destinations it relies on the existing routes G nm , but avoids cyclic pathways that may bring it back closer to s. This resembles pre-planned routes, which, on the one hand, must rely on existing lines, but on the other hand seek the most efficient routing from source to destination, avoiding cyclic pathways. To construct this directed network, we first used the depth first search (DFS) algorithm [9] starting from the therapeutic source node s. At the first step one traverses along all existing G mn links from s to its nearest neighbors, i.e. s's first generation descendants. One then continues to advance along the edges from each of these descendants, seeking additional new second generation descendants, then third generation and so on. In this process the G nm links traversed in each step may be in one of four groups (Fig. S2a): Tree links -linking igeneration nodes to i + 1-generation nodes; Forward links -linking i-generation nodes to i + k-generation nodes (k > 1); Back links -linking i-generation nodes to i − k-generation nodes (k ≥ 1); Cross links -all remaining links, connecting separated branches of the spanning tree. Once we have covered the entire network we eliminate all Back links to obtain the s-dependent directed air-routes network G {s} nm , in which only the routes leading away from s remain. Applying this algorithm to the air-traffic network, starting from s=ITM we arrive at a directed G {s} nm with 19, 605 directed links, out of the original G nm , which comprised 38, 377 links, of which the vast majority were bi-directional, i.e. G nm = G mn .
To obtain the relative fluxes along the directed links of G {s} nm we use the empirical fluxes of F nm to construct the acyclic capturing the volume of (human) travel along each directed route. The available therapeutic at each node m is split among the existing outgoing routes, with flux proportional to F {s} nm . Indeed, the busiest routes, which carry the largest volume of human travelers are also the ones shipping the largest volume of commodities, in this case drug doses. This provides the outgoing fluxes from n to m as allowing us to construct Z nm from (1.41), expressing M n and M m via (1.47).
Balanced distribution (Network 3). In the main text we have shown that the primary lacuna of commodity flow is that, relying on existing air-travel fluxes, it inevitably leads to an unbalanced distribution, in which few nodes receive the therapeutic in excess quantities, while others remain deprived. Hence we consider a balanced network, which calculates for each node, all the downstream destinations along its descending routes, and divides the therapeutics accordingly. This represents an improved, albeit not optimal, physical distribution scheme. Indeed, a fully optimal distribution scheme is extremely difficult to attain, and practically impossible to execute. In fact, even the currently proposed balanced scheme requires a level of control over the air-traffic network that is extremely unlikely to apply in realistic scenarios. Moreover, optimal schemes in general, and our balanced scheme in particular, require us to intervene in the natural flow of commodities, by enforcing external control over the volumes of shipments that travel between destinations, hence they potentially come at a tremendous economic price-tag. Still, we examine this best case scenario, assuming a fully controlled balanced distribution scheme, as a test case for centralized mitigation under a well-designed distribution network.
To construct the balanced network we repeated the process leading to Network 2, constructing G nm from (2.5) and using DFS from the source s. This time, however, we only kept Tree links, leading to an optimized network, in which all target nodes n = 1, . . . , N are linked to the source s via G nm 's shortest paths. Hence we arrive at G {s} nm , in which there is a single directed path Π(s → n) from s to all target nodes n. This path includes a minimal number of flight transfers, as enabled by the existing routes in G nm , hence representing an efficient distribution network, which adheres to the bounds of our current overall aviation capacity.
Each node n is then assigned a branching index K n , quantifying the total population in all of the nodes that are downstream from n. To be specific, consider the path providing the shortest sequence from s to the target node i. For some intermediate node m we define the group (Fig. S2b) which consists of all nodes i whose shortest path from s traverses m, namely nodes that are downstream from m. Note that m itself is also included in K m . Additionally, if m is a leaf node, namely a node with no outgoing links then K m = {m}, including only the node m itself. The corresponding branching index is defined as  Figure S2: Constructing the centralized distribution networks B mn . We constructed three distribution networks. (a) Diffusive (Network 1), comprising all links of the original air traffic network; Directed (Network 2), in which we eliminated Back links (green) to avoid cyclic pathways. Hence shipped doses avoid returning towards the source s (here represented by node number 1). Balanced (Network 3), in which we eliminated also Forward links (red) and Cross links (blue), providing a perfect tree network with a single directed route between s and all target nodes. (b) In Balanced (Network 3) each node n is assigned a branching index K n quantifying the total population downstream from n, i.e. sum of M i over all nodes i whose path from the source 1 traverses n. For example, the nodes downstream from 6 are 9, 10 and 6 itself, hence K 6 = M 6 + M 9 + M 10 = 10. To achieve a balanced spread of the therapeutic, each node must receive doses in quantities proportional to its branching index, to serve all its downstream target nodes. Therefore, in Network 3 the incoming link to n is proportional to K n (Eq. (2.11)).
namely the total population in all nodes located downstream from m. Equation (2.10) captures the amount of doses that must be delivered from s to m if the entire population is to be served. The flux from n to m is then designed to be proportional to K m , namely ensuring that of the available doses q n (t) at n, the share that reaches m is balanced with m's projected downstream shipments, amounting to K m . As before Z nm is constructed from (1.41).
In Fig. S3a we present Res vs. C as obtained for Networks 1 -3 (orange, yellow, brown). For comparison we also show the results of decentralized mitigation (green). As expected decentralized mitigation achieves consistently better results compared to all three centralized distribution schemes. Among the centralized schemes the balanced Network 3 provides the best results, as expected by its egalitarian design. Interestingly Notwork 1, incorporating un-directed diffusive distribution is superior to the directed Network 2. The reason it that in the absence of cycles, a large volume of the available therapeutic is accumulated at the final destination nodes (leafs), becoming effectively lost. The diffusive Network 1 allows the drugs to continue cycling until reaching a destination where they are needed, and hence consumed. The Gini index is presented in Fig. S3b, showing a similar trend among the three centralized distribution strategies.

The rate matrix ξ sn
Consider the spread of a therapeutic under centralized mitigation via Eq. (1.59), namely To track the supply distribution across all nodes let us ignore the local consumption, i.e. first term on r.h.s., and, for simplicity, set t R = 0, and the normalized capacity to κ s C s = 1. Using our diffusive spread of Eq. (2.3) we arrive at in which one unit of the therapeutic is produced at s per unit time, then diffuses through the distribution network B mn to reach all other nodes. To asses the rate of the incoming therapeutic at some destination n consider a path Π(s → n) = {s − − → n}, comprising a sequence of nodes and links along the drug distribution network from s to n. Each link B ij in this pathway represents the fraction of doses present in j that will be routed to i, i.e. B ij = B i←j . Hence a unit dosage exiting s that travels to n through Π(s → n) will be diluted to This represents the rate by which Π(s → n) delivers the therapeutic to n. In case there are multiple paths between s and n, the total rate ξ sn can be approximated by the sum over all routes as For routes of length l this can be expressed as a matrix product ξ {l} sn = [B l ] sn , i.e. the s, n entry of B l . Therefore, the total rate of incoming therapeutics into a target node n is the sum over the contribution of all individual pathways, expressed by as appears in Eq. (6) of the main text. Taking B nm from the appropriate network (Net-works 1 -3) we obtain the rates ξ sn , providing the amount of doses q n reaching destination n per each unit shipped from the source s. The density function P (ξ) captures the fraction of nodes whose rate is ξ sn ∈ (ξ, ξ + dξ). In Fig. S3c -e we present P (ξ) as obtained from Networks 1 -3, setting the source node to be s =ITM. We find that in our air-travel network the longest path has five legs, therefore we set L Max = 5 in (3.6). A node n is considered saved (S n = 1) if, under treatment, its mortality rate is reduced to at least half of that observed in the absence of treatment, namely where r U n (t) is extracted from Eqs. (1.53) -(1.60) and r n (t) from the untreated (1.2) -(1.4). The probability of node n to be saved is strongly dependent on the availability of the therapeutic at n, and hence on its rate ξ sn . To observe this we show in Fig. S3c e the probability P (S ∩ ξ) for a random node n to have ξ sn ∈ (ξ, ξ + dξ) and also have S n = 1. As predicted, we find that the uneven dilution strongly impacts the inequality in the S-probability, where nodes with large ξ sn have a much higher probability of being saved than those with small ξ sn . The effect is less pronounced in Network 3, which is specifically designed to balance the inequality in B nm .

Network pathways
Consider a random symmetric binary network G nm = G mn with an arbitrary degreedistribution P (K), and randomly selected symmetric weights W nm = W mn . Combining these two matrices we construct a model transportation network F nm = G nm W nm , resulting in a weighted network that is extracted from the configuration model ensemble [10]. We can now construct the consequent distribution network B mn , whose weights capture the volume of commodity flow from n to m, via Eq. (2.3) as In (3.8) the weighted degrees S m = N n=1 B mn are normalized to unity, as in Network 1 of Sec. 2.2, capturing the fact that all incoming commodities must also flow (or be consumed) out to other destinations, hence the total out-flux amounts to unity.
We define the un-weighted degree of m as counting the number of routes in/out of m, allowing us to express S m as 3. The effect of heterogeneity is most pronounced in Network 2, consistent with its lowest performance levels among the three networks. Note, that in Network 2 the power-law P (ξ) is truncated for large ξ. As expected, heterogeneity is least pronounced in Network 3, which is designed for balanced distribution.
where B m is the average weight of m's K m outgoing links. Using the fact the S m = 1 we can extract the average weight as hence for a node with K m outgoing routes, we have each route carrying, on average, 1/K m of the available commodities for further shipping. Next we consider m's residual degree capturing the average un-weighted degree of m's nearest neighbors. Averaging over all nodes (m = 1, . . . , N ) we obtain the network's mean residual degree which expresses the average degree over all nearest neighbor nodes. This average neighbor degree is also commonly expressed as [10] where K n represents the n-th moment of P (K). Equation (3.11) allows us to evaluate the average weight of a link associated with nearest neighbor nodes as The analysis above of the residual degree K helps us characterize the pathways in G nm , and hence, through (3.5), obtain the distribution of the rates ξ sn . Indeed, the average node has K nearest neighbors, which, in turn, have, on average K neighbors of their own. Therefore, a node has, on average, K K next neighbors and K K 2 neighbors at distance l = 3, and so on. For a sufficiently large network, with N → ∞, this predicts that the number of nodes at distance l from a randomly selected node grows as N (l) = K K l−1 , an exponential growth that we approximate by [10] where λ = ln K . Equation (3.16) is valid as long as l λ −1 ln N , capturing the network's average path length l , beyond which the shells N (l) cease to grow exponentially. It follows from this equation that starting from an arbitrary source node s, the probability of a randomly selected node n = 1, . . . , N to be at distance l sn = l is a valid approximation as long as l ≤ l .

Rate distribution P (ξ)
In a large random network, due to the local tree-like structure, there is often a single shortest path Π(s → n) between a source s and its target n, whose length is l sn . The typical nodes along such path are, by definition, nearest neighbor nodes, and hence, on average their degree is K , and the l sn weights linking between them are B . Consequently, the rate ξ sn in (3.5) can be approximated by Therefore, on average, we have rates decaying exponentially with network distance l, as where, as above (Eqs. (3.16) and (3.17)) λ = ln K . Equation (3.19) predicts that the average node at distance l from s receives the therapeutic at a rate proportional to e −λl , portraying an exponential dilution as the drugs propagate downstream from s. We can now use (3.19) to obtain the density function P (ξ), by transforming from the random variable ξ to the variable l, whose distribution is given by P (l) in (3.17). Hence we write Taking l from (3.19) we have l = −λ −1 ln ξ and hence | dl/ dξ| ∼ ξ −1 . Using (3.17) and (3.19) we can also write P (l) ∼ ξ −1 , hence, together we arrive at where ν = 2, precisely the prediction of Eq. (7) in 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. However the power-law form of (3.21) is expected under almost any network construction. Indeed, as long as the network exhibits an exponential growth of the form (3.17), and the downstream dilution exhibits an exponential decay of the form (3.19), we expect a scale-free P (ξ). To observe this, consider, for example that the two exponential functions exhibit different growth/decay rates, i.e.

Critical capacity C η
Consider the probability density P (ξ). For a finite system of N nodes it has a minimal rate, ξ min , such that Denoting the normalization constant by Z, we express P (ξ) as P (ξ) ∼ Zξ −ν , and use (3.25) to obtain The maximal rate for a system of N nodes can be obtained from P (ξ ≥ ξ max ) = 1/N , namely, the probability to exceed ξ max is 1/N , providing and in turn At the same time, we can independently evaluate that ξ max = ξ ss , namely the rate at the source itself, which is by definition equal to unity, and therefore For example, under ν = 2, we find that ξ min , which represents the rate at the majority of the nodes is approximately 1/N , an extremely low rate in case N → ∞.
Next we seek the critical capacity C η to save an η-fraction of nodes. Consider a node n with rate ξ sn . For each unit therapeutic produced at s, n receives a fraction ξ sn , and hence under capacity C s , the supply rate at n is of the order of C n = ξ sn C s . The meaning is that the therapeutic availability at n, ignoring local consumption and just focusing on the supply from s, is, roughly, q n (t) ∼ C n t. For n to be saved, i.e. S n = 1, we require C n ∼ 1, such that n is fully supplied (q n (t) = 1) within a finite time. Therefore, the critical ξ sn for saving n is approximately To save an η fraction of all nodes we set C s in (3.31) to C s = C η , namely a critical rate of ξ Crit = 1/C η , and require that all nodes with ξ sn ≥ ξ Crit are saved. Therefore, we write setting the fraction of nodes whose rate exceeds the saving threshold ξ Crit to η. We obtain which, using (3.26), provides where in the last step we used (3.30) to obtain ξ min 's scaling with N . Finally, we use (3.31) to convert ξ Crit to the consequent critical capacity, yielding where φ = 1/(ν − 1), as appears in Eq. (8) of the main text.

Mean supply time T
As explained above, under a given capacity C s , node n, whose rate is ξ sn receives a C n = C s ξ sn fraction of its demand per unit time. Therefore its supply time T Supp,n to fill up its demand (i.e. q n (t) = 1), follows This allows us to transform from P (ξ) to P (T ), capturing the probability density to observe T Supp,n ∈ (T, T + dT ). We write P (T ) dT = P (ξ) dξ dT dT, (3.37) and, using P (ξ) ∼ ξ −ν , ξ ∼ 1/T , we arrive at P (T ) ∼ T ν−2 . As before we seek the normalization factor Z T via which provides The final approximation, relying on T max T min , is relevant as long as ν > 1, a condition which is anyway necessary for P (ξ) to be a valid (normalizable) density function. The mean supply time is given by once again -neglecting the T min limit. Taking Z T from (3.39) we arrive at T ∼ T max , indicating that the mean supply time is governed by the latest supplied nodes. Indeed, given the scale-free nature of P (ξ), the low ξ (and hence large T ) nodes represent the majority of the destinations, and therefore they are the ones that determine the average supply times of the therapeutic. To estimate T max we recall from (3.36) that it the inverse of ξ min . This allows us to use the scaling relationship of (3.30) to write In centralized mitigation each outbreak is characterized by the disease origin at node o and by the therapeutic source from node s, where drugs are produced and shipped globally. In the main text we arbitrarily selected o to be the Bujumbura International Airport (BJM) in Burundi (Fig. S4a, red star), and s to be the Osaka International Airport (ITM) in Japan. In Fig. S4c -k we show Res and the Gini index vs. C as obtained from centralized vs. decentralized distribution for a selection of alternative sources. As expected, the specific source has little impact on the therapeutic efficiency, with decentralized consistently outperforming centralized. Ideally, one expects the best outcome to be observed if s and o are close to each other. Indeed, at the early stages of the outbreak the nodes in the vicinity of o are most severely impacted, and hence if these neighbors of o are efficiently treated, by, e.g., setting s in close proximity, this may lead to a meaningful reduction in global mortality. This goal is best served if s is directly linked to o's neighborhood, in which case the local spread of the therapeutic is most efficient. To observe this we measured Res vs. ξ os , providing an effective distance measure, along the paths from s to o. Setting C = 0.2, we find that, indeed, mortality is higher when ξ os is small (Fig. S4b), however the effect is marginal.

The effect of disease parameters
In the main text we examine the effect of centralized vs. decentralized distribution under the extreme conditions of a highly fatal epidemic, in which R 0 = α/β = 10, conditions that, untreated, lead to an almost complete demise of the human population, i.e. r n (t → ∞) → 1 for all n. In Fig. S5 we examine the results under varying disease parameters, freezing the recovery rate at β = 0.2 day −1 , and setting α such that R 0 = 2, 3, . . . , 10.
As expected, we find that in all cases, decentralized mitigation consistently outperforms the centralized alternative. Note, that the parameter R 0 not only impacts the level of contagion, but also affects the propagation times of the epidemic. This fact touches on the relevant selection of our response time t R , which must be adapted as R 0 is changed. Indeed, setting e.g., t R = 12 days, as we do for R 0 = 10, would represent an unrealistically early response under R 0 = 2, which, over the course of 12 days remains still in its early embryonic stage, likely not even detected as a potential threat. Therefore, to conduct a controlled examination of the impact of R 0 , we update our response time t R as we change the disease parameters. We consistently select t R to be the time t when the global infection level I(t) has reached an η fraction of its peak value, setting η = 0.1 (Fig. S1).

Overview
To examine the potential efficiency of the therapuetic distribution under centralized mitigation we used the commodity flow framework [11] to design an 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, as given by the initial demands d n (0), to all target nodes n = 1, . . . , N . Each link is characterized by a nonnegative carrying capacity C 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 as efficiently as possible, i.e. in minimum time, within the shipping constraints provided by C nm . An important consideration that arises from our earlier discussion is that distribution inequality can severely impact the efficiency of a therapeutic in battling the epidemic. Therefore we wish to distribute the therapeutic to all destinations in parallel, avoiding serial strategies that first fill up selected destinations, then continue shipping to others. To ensure such strategies are avoided we implemented a day-by-day optimization strategy, in which we maximize the commodity flow in each day t = 1, . . . , T . Hence, at each day t we fill up as much of the current demand of all nodes as possible within that single day. We then update the remaining demand d n (t) → d n (t + 1) and reiterate the process at day t + 1, until all demands are satisfied at t = T max = max N n=1 (T Supp,n ).  Figure S5: The effect of disease parameters. We tested decentralized (green) vs. centralized (yellow) mitigation for a range of disease parameters R 0 = 2, . . . , 10, finding that under all conditions decentralized outperforms centralized mitigation.

Maximum flow optimization
At t = 0 we select a source node s, and set the initial demands of all other nodes to d n (0) = M n , i.e. each node's demand is equal to its population size. In each day the therapeutic is distributed from s, aiming for maximum flow [12][13][14], namely seeking to maximize the number of doses exiting s and spreading to all other nodes. Denoting the volume of doses shipped from m to n on day t by F nm (t), this optimization scheme translates to Maximize F s (t) = N m=1 F ms (t) − F sm (t) . (5.1) Using this notation the net supply of doses entering a target node n throughout day t is −F n (t) = N m=1 [F nm (t) − F mn (t)], and hence at the end of each day the demand is reduced by precisely that amount (note that the minus sign in −F n (t) represents incoming flux to n, a counter-intuitive notation, which is motivated by the fact that F s (t) -positive flux -represents drugs exiting the source). Consequently, the n-demand at day t is given by d n (t) = M n + t−1 τ =0 F n (τ ), (5.2) starting at d n (0) = M n , and reducing the demand due to the cumulative therapeutic influx as time advances. The maximization of (5.1) is subject to two constraints that must be satisfied at all t: The first constraint (5.3) ensures that all daily shipments from m to n do not exceed the carrying capacity C nm ; the second constraint (5.4), with d n (t) taken from (5.2), restricts nodes from accumulating therapeutics in excess of their remaining demand, therefore bounding the total incoming doses at day t, i.e. −F n (t), by the current demand d n (t).
To solve the maximum flow optimization problem of (5.1) -(5.4) we used an iterative approach (Fig. 4a, main text): at each time step t we use linear optimization to maximize F s (t) under the two constraints. This provides us with day t's optimized shipping strategy, detailing the daily fluxes F nm (t) through all routes C nm . With these fluxes we calculate the net influx into each node −F n (t), and update the demands d n (t + 1) in (5.2) accordingly. We then use these updated demands to obtain the shipping strategy for the next day t + 1, repeating the process until all demands reach d n (T Supp,n ) = 0 at the final step t = T max . Note that in this process we disregard the flight times themselves, which, if accounted for, might incur additional delays. Instead we assume that the daily time unit accounts also for the lag of each route [15].
Implementation. To construct the carrying capacities C nm we used the passenger fluxes F nm , as obtained from the air travel data [8]. 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. This allows us to write  Figure S6: The effect of the cargo carrying capacity. We constructed C nm in (5.5) using three different values for the passenger/cargo ratio ω, and used the maximum flow optimization of (5.1) -(5.4) to obtain the supply rate distribution P (ξ). As expected the average rate increases as ω is increased, however the fat-tailed (power-law) nature of P (ξ) remains unaffected, indicating that distribution inequality continues to govern the therapeutic dissemination.
C nm = ωF nm , (5.5) where ω is the number of doses that can be shipped as cargo per traveling passenger. Our results are not sensitive to the specific value of ω, as any re-scaling of C nm (capturing doses/day) is equivalent to an arbitrary re-scaling of the time units, having no bearing on the essential findings of Fig. 5 (main text). To select a realistic value, we assume each dose to have a gross weight of ∼ 0.1 kg, which, taking the average passenger mass at ∼ 100 kg, provides a ratio of ω ∼ 10 3 . Therefore we selected ω in the range 10 2 − 10 4 (Fig. S6), with the results in the main text obtained for ω = 10 3 . The initial demands were set to d n (0) = M n , extracting M n from the data via (1.47).