Optimizing hospital distribution across districts to reduce tuberculosis fatalities

The spatial distributions of diverse facilities are often understood in terms of the optimization of the commute distance or the economic profit. Incorporating more general objective functions into such optimization framework may be useful, helping the policy decisions to meet various social and economic demands. As an example, we consider how hospitals should be distributed to minimize the total fatalities of tuberculosis (TB). The empirical data of Korea shows that the fatality rate of TB in a district decreases with the areal density of hospitals, implying their correlation and the possibility of reducing the nationwide fatalities by adjusting the hospital distribution across districts. Approximating the fatality rate by the probability of a patient not to visit a hospital in her/his residential district for the duration period of TB and evaluating the latter probability in the random-walk framework, we obtain the fatality rate as an exponential function of the hospital density with a characteristic constant related to each district’s effective lattice constant estimable empirically. This leads us to the optimal hospital distribution which finds the hospital density in a district to be a logarithmic function of the rescaled patient density. The total fatalities is reduced by 13% with this optimum. The current hospital density deviates from the optimized one in different manners from district to district, which is analyzed in the proposed model framework. The assumptions and limitations of our study are also discussed.

therein increases, which is an important point demanding a quantitative explanation and leads us to expect that relocating hospitals across districts may reduce the total fatalities of TB nationwide.
The optimal distribution of hospitals across districts minimizing the total TB fatalities depends on the concrete form of the fatality rate as a function of the hospital density, which is, however, unknown; The empirically observed negative dependence cannot give this information, as districts are different not only in the hospital density but also in various other properties such as area or population. To address the district-dependent fatality rate, we take a modeling approach, in which the fatality rate is assumed to be identical to the probability of a patient not to visit a hospital and get the medical treatment for the duration period of TB. This is motivated by the expectation that a patient is very likely to be cured once she/he gets a proper treatment in a hospital, given the high success rate of the TB chemotherapy equally applicable to all districts. In this framework, the fatality rate turns out to be an exponentially decaying function of the hospital density, and we are able to derive the optimal hospital densities in all districts the collection of which decreases the total fatalities of TB by 13% from the current value. The predicted optimal hospital density is given by a logarithmic function of the rescaled patient density. Our results delineate an analytic approach to the facility optimization problem under an objective function from a public health perspective. The limitation and further generalization of the results will be discussed.

Results
TB fatality rate and hospital density: Empirical data. The incidence and mortality of TB are well recorded in Korea. In Statistics Korea 20 , we obtain for district = … = i I 1, 2, , 228 in year 2014 the number of the newly reported TB patients N i , the number of dead TB patients (fatalities) D i , the number of private general hospitals H i , and the area a i . Here "district" includes three distinct units for administrative division, Gu, Gun, and Si, with the population ranging from 10 4 to 10 6 and smaller than the metropolitan cities.
We are interested in the fatality rate φ i of TB, defined as the ratio of the number of dead TB patients to the number of new TB patients reported for one year in each district i, i i i It is quite different from district to district, ranging between . 0 01 and . 0 25, as shown graphically in Fig. 1(a). What drives such difference in the TB fatality rate? Taking regularly medical treatments and examinations in hospitals may be the most important for curing TB, which is available in the easy-to-frequently-access medical environment established in the local community. Therefore a difference in the abundance and accessibility of Figure 1. Distribution and relations of the TB fatality rate φ, the areal density of hospitals η, and the patient density ρ in Korean districts at the level of Gu, Gun, and Si. The unit of η and ρ is km −2 . (a) The TB fatality rate φ is represented by color in 228 districts. Seoul, the capital city, has 25 Gu's and is shown separately. (b) Fatality rate φ versus hospital density η. Open circles are the raw data for all districts and filled squares represent the average fatality rate for each given hospital density with the standard deviations as errorbars. The Pearson correlation coefficient is −0.26 with P-value 0.000070 for all districts and −0.25 with = .
P 0 0025 for the districts with both η and φ non-zero. (c) Hospital density η versus patient density ρ. As in (b), open circles and filled squares represent the raw data and the average value, respectively. The solid line fits the averaged data and its slope is . ± .
www.nature.com/scientificreports www.nature.com/scientificreports/ hospitals in the patients' residential districts will be a major factor giving rise to such variation of the fatality rate with district. In this light, we investigate the relation between the fatality rate φ i and the areal hospital density i i i in unit of km −2 . In Fig. 1(b) φ i tends to decrease with η i ; The larger the hospital density is, the smaller the fatality rate is. This correlation is significant with < − P 10 4 . Yet the dependence does not look so strong as expected. This will be shown to be due to that the fatality rate of a district may depend not only on the hospital density but also on other characteristics. In this work, all correlation values are measured in linear scales.
What principle underlies the current spatial distribution of hospitals? The scaling behavior with respect to the patient density has hinted at the answer 14 . The hospital density scales with the areal TB patient density ρ ≡ i i in which α = . ± . 1 05 0 07 when all districts are included, and α = . ± . 0 80 0 045 when the districts having no private general hospital are excluded [ Fig. 1(c)]. Many other properties also scale with respect to the patient density. The patient density is almost linearly related to the population density ρ′ = P A / i i i with P i the number of people living in district i [Fig. S1]. The exponent α for the hospitals in United States is close to 1, rather than 2/3 14 . These results suggest that the profit maximization affects the hospital distribution. For self-containment, let us sketch the corresponding optimization calculations. The sum of the economic profits of all hospitals distributed across I′ districts in a country is given by with ω x ( ) the expected profit of a single hospital having x patients available. On the other hand, the sum of the social costs, such as the travel distances, of patients is given by 1/2 being the expected travel distance of a patient residing in a district of x area per hospital. Then, for a fixed total number of hospitals Our question is then whether the current hospital distribution, seemingly maximizing the economic profit, is the best also for minimizing the total fatalities of TB Can E fatalities be reduced by the redistribution of hospitals across district, i.e., some change of η { } i ? To answer this, we should formulate the total fatalities in Eq. (7) as the objective function and minimize it with respect to the hospital density for the given total number of hospitals in Eq. (6). The fatality rate φ i should be some function of the hospital density η i . If the optimal hospital densities η { } i (opt) are obtained by this optimization computation, we will be able to evaluate the quality of the current spatial distribution of hospitals regarding its capacity of TB treatment. Also we will see immediately how to redistribute the hospitals to reduce the TB fatalities. In the present study we do not consider a variation in the numbers of TB patients N { } i but take them for given; The onset and spreading of the TB or a general epidemic disease depend strongly on the topology of human contact networks and the infection rate, which is another important research topic and has been studied extensively 21,22 . Fatality rate as a function of hospital density: Model. The empirical fatality rate φ i in Eq. (1) can be considered as the probability of a TB patient to die, losing the opportunity to get proper medical treatment in time. Our idea is to approximate the latter by the probability that a patient does not visit any hospital in her/his residential district for a given period = t 3 (TB) years, the empirically reported period of TB duration from onset to either cure or death 23 . In this model framework, it determines the fate of a TB patient whether she/he visits a hospital or not for the period of t (TB) . The patient will recover if yes, but will be dead otherwise. One can see that this is a trapping problem 24 from the viewpoint of a patient; Once a patient (walker) reaches a hospital (trap), she loses the status of a patient (absorbed at the trap). The probability of a walker to survive during a given number of steps corresponding to t (TB) in this trapping problem is translated into the fatality rate of a TB patient in reality. (2020) 10:8603 | https://doi.org/10.1038/s41598-020-65337-x www.nature.com/scientificreports www.nature.com/scientificreports/ Suppose that H traps are uniformly and independently distributed in a two-dimensional Euclidean lattice of × L L sites and that a walker walks around the region, who disappears on reaching any one of the traps. Then the probability of the walker to survive (not to reach any of the traps) after τ steps is given by is the density of traps and τ S( ) is the number of distinct sites visited up to τ steps.  represents the average over different realizations of walks. In the limit λ  reachable when the trap density is sufficiently low or the number of steps is small enough, the survival probability φ can be approximated in terms of the first cumulant of the probability distribution of S as 25 S ( ) which is the exponential function of the trap density λ. It seems that Eq. (9) allows us to relate the hospital density and the fatality rate. However the dimensionless quantities λ and τ S ( ) are not directly available. In random walks in two dimensions, the expected number of distinct visited sites is known to be 24 which is inserted into Eq. (9) to give log (11) with the coefficient = . . c 3 5/1 13 2 known numerically 26 . In the opposite limit of the walker is governed by the probability of a large trap-free region to be formed, which leads to a stretched exponential form φ λτ ∼ log [26][27][28] .
The exponential decay of φ with λ in Eq. (9) holds when the hospital density is sufficiently low. The randomness of the mobility pattern is assumed in obtaining Eq. (11). We should remark that the human mobility pattern revealed by tracing the travel routes of bank notes 29 or the mobile phone records 30 displays deviation from random walk; The radius of gyration of individual trajectories grows logarithmically with time 30 , in contrast to the square-root scaling in the conventional random walk, and such slow diffusion is known to arise under the memory effect [31][32][33] or the spatial quenched disorder 24 . The assumption we make about the human mobility pattern is that the coarse-grained trajectories of individuals on the time scale of = t 3 (TB) years, much longer than the previous studies, show the survival probability given in Eq. (11) like random walks. The coarse-grained trajectory is obtained by neglecting the spots swiftly passed by and connecting the remaining notable places which an individual visits and stays for a while in, such as her/his house, workplace, parks, stores, banks, oil stations, and hospitals. In our model, we are interested in whether a hospital is included in the list of such notable places. We cannot check directly the validity of Eqs. (9) and (11), however, we will present indirect evidence that they are reasonable assumptions.

Lattice constant and dimensionless quantities.
To relate the survival probability in the 2D trapping problem to the fatality rate of TB, we need to convert the empirical data into dimensionless ones of Eq. (11). To this end, we discretize the region of each district i by introducing the lattice constant a i , corresponding to the typical length of one single step or the average distance between adjacent notable places appearing in the coarse-grained trajectories. Then the district is represented by the × L L years. Then the number of steps taken in her/his coarse-grained trajectory for t (TB) in a district i will be given by Plugging Eqs. (12) and (13) into Eq. (11), we find the fatality rate represented as Assuming the validity of Eq. (14), one can estimate the characteristic hospital density η  i by using the empirical data of the fatality rate φ i and the hospital density η i in Eq. (14) as (14) with the estimated coefficient η  i for selected districts are shown in Fig. 2(a). η  i is different from district to district, growing with the patient density [ Fig. 2(b)], which underlies the weaker decay of the fatality rate with the hospital density [ Fig. 1(b)] than would be expected if η  i were identical for all districts. The estimated η  i is the characteristic constant of each district and will be used throughout the optimization computation.
The lattice constant can be obtained by using the estimated η  i in Eq. (15) and solving for a i . Let us denote the solution by  a ( ) i (F) (TB) . To validate it, we compare it with another estimate independent of the empirical values of the fatality rate or the hospital density. We use the data of the number of business buildings B i in each district 20 . The business buildings, including hospitals, are the candidates for the notable places included in the coarse-grained trajectories. The typical distance between adjacent business buildings can be a candidate for the lattice constant, which is given by  Fig. 3(a)]. It corresponds to the annual traveling distance 3300 km which is reasonably close to the empirical value 8478 km of Korea 34 . In Fig. 3(a), the two lattice constants  = ⁎ a a ( ) (F) ( F) (TB) and a (B) show good agreement in their magnitudes, supporting the validity of the assumptions of our model and its formulas, Eqs. (14) and (15). Due to this agreement and Eq. (17), we can see that a large or small value of a i (F) originates from the sparse or dense business buildings in district i. With a i (F) , the dimensionless hospital density λ i and the number of steps τ i taken for t (TB) can be evaluated by Eqs. (12) and (13), which are plotted versus the patient population density in Fig. 3(b,c), respectively. In contrast to the real hospital density η i , λ i is lower in a district with higher patient density [ Fig. 3(b)]. It is attributed to the smaller lattice constants in the districts of higher patient densities, arising from the denser buildings. On the other hand, the number of steps taken for t (TB) increases as the patient density increases, increasing the chance to visit hospitals. To sum up, effectively less hospitals are distributed but the patients take more steps for the given period t (TB) in the higher-populated districts, which explains the slightly lower fatality rates therein than in lower-populated districts as shown in Fig. 1(b).
Optimal hospital density. The fatality rate formula in Eq. (14) is applicable to I s districts having non-zero η and φ in the empirical data. Then one can minimize the total fatalities in those I s districts with z being the Lagrange multiplier. Consequently the optimal hospital density is found to be and the optimal fatality rate is The Lagrange multiplier z is computed by inserting Eq. (19) into Eq.
. We call the hospital density in Eq. (19) optimal in a sense that it minimizes the objective function given in Eq. (18). If the objective function is changed, the optimal density may be changed, which is discussed in the Summary and Discussion section and in the Supplementary Information (SI). Equations (19) and (20) are the main results of the present study. Remarkably the optimal hospital density and the patient density are rescaled commonly by η  i and then related to each other logarithmically. In Fig. 4(a), the arrangement of the data points of the optimal hospital densities on a straight line is contrasted with the scattered distribution of the current (empirical) hospital densities in the (ρ η η η   / , / ) plane in semi-logarithmic scale. The same phenomenon is observed for the fatality rate; the empirical fatality rates φ i 's are scattered but the optimized fatality rates lie on a straight line in the (ρ η φ  / , ) plane in logarithmic scale as shown in Fig. 4(b). The rescaled patient density ranges between 30.80 (Yeonggwang-gun) and 2646 (Songpa-gu), and is larger than .  z 12 9 and thus guarantees η for all i in Eq. (19). More plots of the optimized hospital densities and fatality rates are given in Fig. S2.
The scattered distributions of the empirical data in Fig. 4 show clearly the deviation of the current distribution of hospitals from the optimum minimizing the total fatalities of TB. The minimized total fatalities E fatalities (min) obtained from the optimal hospital distribution is  www.nature.com/scientificreports www.nature.com/scientificreports/ which is smaller than the current value, 1718, by 13%. For this optimization, hospitals among a total of = H 328 total are relocated. In the zero-temperature Monte Carlo (MC) simulation in which some small amount ∆H of hospitals are moved between randomly selected districts only when the attempted relocation reduces E fatalities , the theoretical predictions in Eqs. (19) and (21) are realized in the steady state as long as ∆H is sufficiently small [ Fig. 5(a)]. If H i 's are restricted to be integers (∆ = H 1), the total fatalities in the steady state is 1599.45. The stationary-state results remain unchanged in the simulations with different initial configurations or with gradually cooling down the temperature. For more details of the simulations, see Methods.
To achieve such reduction in the total fatalities, the hospital density should be increased in some districts and decreased in others [ Fig. 5(b)]. For instance, Gyeongju-si should have its hospital density 1.61 times larger than the current hospital density but Yeonggwang-gun 0.534 times larger than the current one. In Fig. 6(a), 143 districts are colored blue (red) if the optimal hospital density is larger (smaller) than the current hospital density. Interestingly, the ratio η η (opt) of the optimal to current hospital density turns out to depend strongly on the rescaled patient density ρ η  [ Fig. 6(b)]. It implies that if the rescaled patient density is large in a district and small in another district, it is recommended to move some hospitals from the latter to the former district. The high and low values of η η (opt) of Gyeongju-si and Yeonggwang-gun can be understood in this line, as they have quite different values of ρ η  , 662 and 30.8, respectively. Such a significant correlation is absent between η η (opt) and the raw patient density ρ [ Fig. S3]. The change of the fatality rate shows the opposite trend to that of the hospital density; The districts with large (small) rescaled patient density ρ η  have their fatality rate decreased (increased) as they gain (lose) hospitals [ Fig. 6(c)]. It is remarkable that the changes of the hospital density and the fatality rate are determined not by the patient density or the population density but the rescaled patient density.

Summary and Discussion
We have here proposed a modeling framework which predicts the optimal distribution under a general objective function, going beyond the previous descriptive explanations for the facility distribution. In deriving the optimal hospital distribution over districts for minimizing the TB fatalities in the whole country, we have found that the characteristic hospital density of each district plays an important role in the optimization. The random-walk nature of the coarse-grained trajectories of patients has been assumed in establishing a theoretical model, and the lattice constant of each district has been introduced to connect the theoretical results and the empirical data. The incorporation of such heterogeneity of districts in the theoretical study of the facility optimization is done only in the present study and can be useful in future studies.
Examining the assumptions and limitations of the proposed model may help better understand and improve its predictive power. The exponential decay of the fatality rate with the hospital density given in Eq. (14) is valid when the dimensionless hospital density is low, λτ τ  1 log 26 . The empirical data analyzed in the present study stay in this regime; λτ τ log ranges between 0.23 and 0.37. If we were to extend to the case of high hospital density or the hospital locations being no more independent of one another, the fatality rate might behave differently from Eq. (14).
( ) a decreasing convex function such as exponential, stretched exponential, or power law, the optimal hospital density will be given by = − ′ x . The case of f x ( ) being a power law is presented in the SI. Regarding the robustness of the functional form of the fatality rate, it will be of interest to investigate which model of walk with traps exhibits such a power-law survival probability.
We have counted only the private general hospitals, but there are mostly found one or two public health centers in a district, which can also provide the medical treatment to TB patients although its portion may not be large 35 . One can optimize the private hospital distribution considering the contribution of the public health centers to the TB treatment, which is presented in the SI and Fig. S4. The results remain unchanged qualitatively. Extending the trajectories of patients to the nearby districts can be one way of making the model more realistic, which will address how the similarity or dissimilarity of adjacent districts may affect the fatality rate and the total fatalities. A correlation may be expected between the economic level of a district and the fatality rate of a disease, which is not significant in case of TB in Korea [ Fig. S5], as the TB treatment is covered by the national medical insurance in Korea 36 , but should be considered in the application to other diseases or other countries. When a given number of hospitals can be opened or should be shut down for financial or other reasons, our results will be helpful for the investigation of the optimal locations. . Upper and lower triangles are used for the data points with the optimal hospital density larger and smaller, respectively, than the current one. The filled square is the average and the errorbar is the standard deviation of the ratio η η (opt) . (c) Plot of the ratio of the optimal to current fatality rate φ φ (opt) versus ρ η  .