Quantification of the resilience of primary care networks by stress testing the health care system.

There are practically no quantitative tools for understanding how much stress a health care system can absorb before it loses its ability to provide care. We propose to measure the resilience of health care systems with respect to changes in the density of primary care providers. We develop a computational model on a 1-to-1 scale for a countrywide primary care sector based on patient-sharing networks. Nodes represent all primary care providers in a country; links indicate patient flows between them. The removal of providers could cause a cascade of patient displacements, as patients have to find alternative providers. The model is calibrated with nationwide data from Austria that includes almost all primary care contacts over 2 y. We assign 2 properties to every provider: the "CareRank" measures the average number of displacements caused by a provider's removal (systemic risk) as well as the fraction of patients a provider can absorb when others default (systemic benefit). Below a critical number of providers, large-scale cascades of patient displacements occur, and no more providers can be found in a given region. We quantify regional resilience as the maximum fraction of providers that can be removed before cascading events prevent coverage for all patients within a district. We find considerable regional heterogeneity in the critical transition point from resilient to nonresilient behavior. We demonstrate that health care resilience cannot be quantified by physician density alone but must take into account how networked systems respond and restructure in response to shocks. The approach can identify systemically relevant providers.

For the last 50 y, health-related expenditures in almost all western countries have been growing faster than national incomes (gross domestic product) (1). This has raised concerns about the sustainability of health care systems all across the Organization for Economic Cooperation and Development (2). In several developed countries, health care demand will further increase, because the population is aging and the prevalence of chronic disorders is increasing (3). The situation is exacerbated by impending retirement waves (4). Is there a point beyond which these pressures will severely impair the quality of care? If so, how close are we to it? To answer these questions, a quantitative understanding of the resilience of health care systems is required. Resilience quantifies the rate of recovery and the extent to which a system is able to recover from disruptive events (5). In health care systems, such events include sudden increases in patient numbers or reductions in the number of health care providers within a specific region. Resilience captures how fast and the extent to which it is possible to deliver adequate health services to the entire population in the wake of such a shock.

Health care is undergoing a digital revolution (6) driven by the increasing availability of observational health care data (7). As countries adopt systems of shared electronic health records, such data become available at national scales (8). These systems enable analysts to answer questions, such as “Who did what, when, for whom, where and at what costs?,” for practically all medical services in a given country. For instance, in Austria, it has been shown that such data can be used to identify genetic, environmental, and epigenetic disease risks (9, 10) or to investigate how individual health care providers coordinate with each other in the treatment of patients (11–14). Health care providers are embedded in multiple formal and informal relationships, because they share information or patients. These relationships can be formalized in so-called patient-sharing networks, which consist of nodes (providers) connected by links if they share the same patients (15, 16). Such networks show large structural variations that can be related to differences in the cost and quality of care (17–19).

Here, we show how to quantify the resilience of health care systems with respect to changes in local densities of health care providers by means of dynamical simulations of structural changes in the patient-sharing networks. Further, we show how this method can be used to benchmark and stress test 121 Austrian regions (political districts) in terms of their resilience. The central idea of our approach is as follows.

Consider 4 physicians , , , and who share patients with each other: say 20% of ’s patients have also visited , and 10% have seen (Fig. 1). Links between physicians may arise for a multitude of reasons (e.g., because is ’s holiday locum, is on maternity or sick leave, etc.). How doctors share patients is given by a network, , in which doctors are nodes connected by patient-sharing relations. Assume that, at time , doctor is closed for business. As 20% of ’s patients already have a treatment relationship with , it is natural to assume that many of these patients will now seek treatment with . The removal of node induces a displacement flow of patients along the link from to but also from to . By receiving new patients from , both and will get closer to their maximum capacity. In Fig. 1, this is shown by the change in node color. Green means high spare capacity; red means that the capacity limit is reached. In the example, now has exceeded its capacity (received more patients from than can be treated within reasonable time). Doctor must, therefore, in the next time step send the excess number of newly inherited patients to yet other doctors (along the links in the patient-sharing network), here to and (Fig. 1). Nodes and get closer to their limits but are still capable of absorbing more patients. The removal of doctor leads to a cascade of patient displacements of size 2. In other cases, where doctors are closer to their limits, cascades can become large and eventually span a large region of the patient-sharing network.

Fig. 1.

Schematic representation of patient displacement dynamics. (A) Doctors are represented as nodes (size represents the number of patients treated per year). They are linked if they share patients in the patient-sharing network, (black arrows). The color represents their current capacity; green means that they have capacity, and red means that they can no longer accept new patients. (B) Doctor retires at time step 1; his/her patients are distributed to other doctors according to the weights of the links from to and from to (yellow arrows). This, in turn, changes the capacity of the other doctors. (C) As has reached its capacity limit (red), he/she must send patients to other doctors (blue arrows from to and ). This creates a cascade of patient displacements of size 2. D–F show the same steps as in A–C in a simulation of a realistic environment. Doctors are localized (due to data protection) at random locations within a district, and a real patient-sharing network is used. (E) A doctor is removed, and his/her patients are shared (yellow). (F) Those doctors who reach their capacity send excess patients to others in a second round (blue). At this point, all patients are cared for, and the model dynamic terminates.

A highly resilient health care system should be able to redistribute ’s patients with a minimal number of patient displacements in a short interval of time. The initial shock (’s removal) is then quickly absorbed, and the system becomes fully functional again soon afterward (all patients find a new doctor). A nonresilient system, however, is characterized by cascades of patient displacements that push multiple doctors beyond their capacities. If a substantial number of patients do not find a new doctor, the health care system will essentially lose its ability to deliver adequate care. We can identify 2 related indicators to distinguish resilient from nonresilient behavior. The higher the resilience of a health care system, 1) the lower the number of displacements that the removal of doctors causes, and 2) the lower the number of patients unable to find a new doctor. The systemically beneficial doctors (i.e., those who contribute most to regional resilience) are those who take over the largest shares of patients.

Cascading processes are examples of dynamical phenomena that take place on networks (20). To model such processes, a localized perturbation is considered by shocking or removing a single node. This initial event spreads via the links of the perturbed node to other nodes, which might trigger another step in a cascade as those nodes propagate the shock to their neighbors and so on. Such processes can be formulated by means of recursive centrality measures (e.g., the PageRank algorithm) (21) or models that consider load distribution on networks (22). In concrete applications, these network measures often require modifications that reflect specific properties of the system under consideration, such as the propagation of shocks between financial institutions (23), failures in power grids (24), or cascading failures in interconnected infrastructures (25).

Here, we develop a data-driven computational framework to estimate the impact of doctor removals through cascading processes of displacements on patient-sharing networks of practically all physicians in Austria. We construct patient-sharing networks of primary care providers (PCPs) for 121 districts from an extensive dataset containing about 97% of all outpatient contacts over 2 y in Austria (9–12) (). We formulate a dynamical model that simulates the removal of one or several providers and computes the size and duration of the resulting cascades of patient displacements in the following way.

Every PCP is a node in the patient-sharing network with weighted directed links from to . The link strength, , corresponds to the number of patients of who occasionally also visit PCP . In every quarter of a year , every PCP treats unique patients. The average number of unique patients who a doctor sees in a quarter is , where is the total timespan of the data. Every PCP is further characterized by a fixed capacity , which is estimated from historical data. In the simplest case, we assume , where is a free model parameter.

The model dynamic takes place on a timescale, , that is shorter than a quarter, say days. Initially, each patient is assigned to the PCP who he/she most frequently consulted in the past. A PCP is in 1 of 3 internal states: available, fully booked, or unavailable (removed) (). Assume that a PCP is removed from the network of district at time . Those patients who usually visited PCP now transfer to with probability . We allow for the possibility that not all patient displacements follow the links of the patient-sharing network. With probability , patients select a random doctor in the same district with a uniform probability (). To every PCP , we assign the average number of displacements, , that ’s patients must undergo before finding a new and available PCP. Ranking PCPs according to their value of (from high to low) identifies physicians with the largest contributions to systemic risk; we call this rank the CareRank of a PCP. For each PCP , we also measure average systemic benefit, , which is the fraction of displaced patients who end up at averaged over removals of all other providers in the district. The displacements, , and benefits, , are proxies of the systemic risk and benefit contributions of every doctor ; a definition is in .

We use this model to quantify the resilience of individual regions in which multiple PCPs are removed. The set of doctors removed at time is denoted by . The size of this initial shock is the fraction of PCPs who become unavailable at , , where is the number of doctors in district . Following this shock, we count the number of patients in district unable to find an available PCP, (). We refer to as the number of “lost patients.” For each district, , we are interested in the smallest shock size, , for which holds in a certain fraction of model runs. This means that there will be patients no longer able to find primary health care services within a given district. At this critical shock size, the district has reached its “resilience point.” The parameter, , serves as our resilience indicator. We show that, surprisingly, the critical doctor removal density is practically uncorrelated with regional physician densities, a conventional indicator to assess health care coverage (26). To explore how the resilience indicators depend on properties of the PCPs and the networks that they are embedded in, we use 2 different types of linear regression model ().

We consider 4 alternative model variants. First, we assume that doctor capacity, , can be estimated from the historically observed fluctuations in a doctor’s patient numbers (i.e., is proportional to the variance of ). Second, is implemented as a multistep function to take differences in staffing into account (that is, physicians hire additional staff, which increases their capacity by a constant factor). Third, the next variant is equal to the main variant except that patients seeking a new PCP do not contact the same doctor twice during their search (they perform a self-avoiding random walk on the patient-sharing network). Fourth, the dynamics of the main model variant is studied on the countrywide patient-sharing network without being broken down into districts. has a description of these model variants.

Results

The model dynamic is illustrated in Fig. 1. Initially (Fig. 1), all doctors operate well below their capacity (green). At time , PCP becomes unavailable (Fig. 1). His/her former patients seek a new doctor on the patient-sharing network (yellow in Fig. 1). PCP now is fully booked. At , patients can no longer be accommodated by and move from PCP to nodes and (Fig. 1). After all patients find a new PCP, the dynamic stops. A web-based interactive visualization of a simplified version of the model on a real regional primary care network is available online (https://csh.ac.at/vis/med_public/pcn-resilience) (). Structural properties of these patient-sharing networks have been reported previously (11–14).

We now study the validity of 2 central model assumptions. First, we test whether patients who lose their PCP are indeed displaced along links in the patient-sharing network. We identify as removed doctors those who had at least 100 quarterly patients on average in the first year but no patients in the second half of the second year; 28,795 patients with at least 2 different physicians were displaced this way. Of those, 84% most frequently consulted a PCP in 2007 who they had already seen in 2006. In these cases, the removal of a doctor did indeed lead to displacements along the patient-sharing network. Second, we inquire to what extent the nationwide patient-sharing network can be decomposed into individual districts; 77% of links between doctors from the same district are nonzero compared with 3.6% of links between doctors of different districts being nonzero (). How these interdistrict links influence the model results is investigated in the model variant that uses the countrywide patient-sharing network.

Systemic Relevance of PCPs.

We next determine the average number of patient displacements, , caused by the removal of doctor , . As the model is not deterministic, we performed 500 model runs. The median of , the average quarterly patient number per PCP, is 945 (corresponding to about 10 patients per day) (Fig. 2). Fig. 2 shows the relation between and average displacements . The most systemically relevant PCPs cause almost 3 displacements per patient on average, while many cause slightly more displacements than the theoretical minimum of 1. Doctors with displacements close to this minimum tend to have low patient numbers within the range from 20 to 500. The majority of physicians have patient numbers between 500 and 2,500, for which we observe displacements that vary between 1 and 3. These 2 “modes” of the bivariate distribution of physicians in the plane result in a weak linear correlation (Pearson’s , value of ). To explore how differences in relate to other network or demographic properties, we perform univariate and multivariate regression analyses (). Overall, high-impact doctors tend to have high numbers of patients, low numbers of links with high weights, low numbers of closed triadic relationships (low clustering), and low centrality in the network. In Fig. 2, we show that the systemic risk contributions of doctors, , can substantially deviate from their systemic benefit, (Pearson’s , value of ). PCPs in the upper left region of Fig. 2 combine high systemic risk contributions with low benefit, whereas the lower right region shows physicians with high benefit and relatively small risk contributions.

Fig. 2.

Systemic risk profile of Austrian health care providers. (A) The distribution of average quarterly patient numbers of doctors has a median of 945 patients. (B) Displacements, , for every doctor in Austria tend to correlate with patient numbers of doctors (Pearson’s , ). The color encodes the number of doctors with a given (, ) pair. (C) Systemic risk contributions of doctors, , show only little correlation with their systemic benefit (Pearson’s , ), . The 4 quadrants indicate regions where and are above or below their population medians, respectively.

Resilience of Districts.

After removal of a single doctor, patients typically find a new doctor in district , ; no patients are “lost.” In the situation where a larger fraction of PCPs is removed, this can change. We now ask at which critical fraction we find the onset of lost patients, . indicates the location of a regime shift in the model behavior (in the Introduction).

Fig. 3 shows the number of lost patients, , as a function of the shock size for the district of Reutte in Tyrol. Doctors are removed sequentially. We show 2 different sequences (green and blue in Fig. 3). The smallest value of for which becomes nonzero depends on the sequence order. A critical can be defined using the sequence that leads to the largest (upper bound, red arrow in Fig. 3) or smallest (lower bound) (). We perform 500 model runs (50 different choices of specific sequences for 10 model realizations) in which doctors have been removed initially. In Fig. 4, the upper bound for the resilience indicator for each district is encoded in the district color from green (most resilient) to red (least resilient). For most districts, the transition occurs after about of the doctors are removed (). However, there are also districts where the transition occurs for substantially smaller (about 20%) or larger (about 40%) values of . shows that the width of this transition varies substantially across districts. Note that depends on the choice of the capacity parameter and is, therefore, not in itself informative unless reasonable choices are made. However, the relative ranking of individual districts by their for regional comparisons can be carried out for any .

Fig. 3.

Number of patients, , who cannot be cared for as a function of the fraction of unavailable PCPs, , in district . is shown for 2 different scenarios where doctors are removed in a different order: sequence (green) and sequence (blue). Labels display indexes of the removed PCPs at each step (i.e., the green sequence first removes PCP 2 followed by 13, 10, and so forth, whereas the blue sequence first removes PCP 4 followed by 7, 3, etc.). The shaded area envelops all observed values of (100 CI). Sequence gives a scenario where 44% (8pcps) have to be removed before losing patients, whereas 22% (4pcps) are sufficient to put the district in a condition where it cannot care for all patients for sequence . The red arrow marks the position of the critical fraction , which is the smallest such that holds for each observed sequence.

Fig. 4.

Map of Austria that shows the upper bound of the resilience indicator, , for all districts. Districts colored in green (red) have a particularly high (low) resilience: that is, critical removal fractions of .

In Fig. 5, we compare the lower bounds of the resilience scores, , with the de facto standard indicator for health coverage (i.e., physician density; number of PCPs per thousand population). shows a similar comparison using the upper bound of . In both comparisons, districts with similar resilience scores, , can have physician densities that vary by up to 1 order of magnitude. The regression analysis additionally shows a negative correlation of the resilience scores with the district-averaged clustering coefficient [ Pearson’s , ] and a positive correlation with district-averaged closeness centrality, [ Pearson’s , ] (). Both of these correlations are confounded by the demographic properties of the districts ().

Fig. 5.

Resilience vs. (primary care) provider density. Every circle (size is proportional to the district population) is a district with its provider density (number of PCPs per thousand population) on the axis and the lower bound of resilience indicator, , on the axis. While there is some correlation between them (Pearson’s , ; Spearman’s , ), physician density can vary by up to 1 order of magnitude for districts of similar resilience.

Robustness.

We obtain qualitatively similar results in all 4 different model variants and for a wide range of choices in the model parameters (). Results for the patient displacements, , present no qualitative change with respect to the standard variant for values of in the range from 0.01 to 0.1. For even larger values of , cascade sizes approach 1 for many patient displacements, whereas for smaller values, the cascades might easily span the entire system, even for small shocks. We study 2 alternative definitions of the doctor capacity, namely by inferring from the observed variance of patient numbers, , and by assuming a multistep function of capacity to take different levels of staffing into account. The model was also evaluated on the countrywide patient-sharing network (patients can cross districts) and by assuming that patients perform a self-avoiding walk on the network. Due to computational costs, particularly for the latter 2 variants, the results of these variants are compared for the doctor displacements, . Overall, we find substantial correlations between all model variants, in many cases with correlation coefficients close to 1 (). Considering pairwise comparisons between the main model and the other variants (), we find the lowest agreement with the variance definition of doctor capacity (for very low values of ) and with the multistep variant, where we observe correlations around 50%. All other correlation coefficients fall in the range between 70 and 95%. Finally, we confirmed that the relations between our doctor- and district-level systemic risk measures, and , show similar correlations with other demographic and network properties as in the main model ().

Discussion

The primary aim of this paper was to quantify the resilience of regional primary care networks on a fully data-driven basis. We were able to quantify resilience on 2 scales. First, we determined the systemic relevance of individual doctors by estimating the number of patient displacements, , caused by her/his removal. Second, we were able to estimate the critical fraction, , of PCPs who can be removed before the regional health care service breaks down. We developed a full-scale simulation model for how regional patient flow networks reorganize after the removal of one or several doctors. By full scale, we mean that the actual data of each Austrian patient and each PCP are represented as a fully anonymized individual avatar in the model. Avatars are used to infer patient-sharing networks and capacities of doctors. The decisions and behavior of these avatars (that is, how they choose their doctors based on the patient-sharing networks in which they are embedded) have been formulated and calibrated on a large-scale database of observational health care data. The model has 2 relevant free parameters: the shortcut probability, , that captures whether patients choose new doctors through the network, or if they base their choice on other factors, and the capacity parameter, , that quantifies the willingness of doctors to accept additional patients. The introduction of avoids dynamical traps in the model and has only a marginal impact on the results. We showed that our main results are robust with respect to the remaining free parameter, the capacity parameter, . Therefore, it is unlikely that our results are idiosyncratic features of particular choices of the used free parameters, but rather, they reflect genuine structural and dynamical properties of primary care networks.

Regions show considerable differences in their resilience scores. Districts with similar resilience can have very different physician densities (a difference of up to 1 order of magnitude). Physician density assumes that doctors and patients circulate freely in their regions and meet with a probability that is independent of their actual position. This view neglects the actual structure of the underlying patient-sharing networks. Our results show that resilience of the primary care sector is largely a network effect that is determined by the onset of cascading processes of patient displacements triggered by the removal of one or several nodes. After we take the demographic characteristics of the districts into account (e.g., their population, number of doctors, population density, etc.), we find no significant correlations between the regional resilience scores, , and a number of conventional network measures (connectivity, clustering, centrality, and so on). Our measure of resilience, therefore, quantifies a genuine network capacity effect of how efficiently the network distributes cascades of patient displacements. Consequently, policy makers should exercise caution when using physician density indicators to estimate the impact of changes in the density of care providers on health service accessibility or coverage, since ignoring the network structure of the health care system might severely under- or overestimate these impacts.

Alongside the systemic risk contributions of individual doctors, we show how to quantify their systemic benefit, , in terms of how many patients they typically absorb in a patient displacement cascade. We find a large number of PCPs who combine relatively high systemic risk with rather low systemic benefits or low risk with high benefit. This is to be expected, since the first is determined by the flow of patients from the PCP, while the second is determined by the flow to the PCP. As the network is not symmetric, these need not be the same. Our results suggest that the health care system could be made more resilient by protecting doctors with high benefit, (or prioritizing their immediate replacement after they leave). Our study also highlights the necessity for an investigation on how specific structural changes (e.g., increased number of multiprofessional primary care centers) impact resilience.

Several limitations originate from the quality of the underlying dataset; a thorough discussion is found in refs. 11 and 12. For resilience assessment, a relevant limitation is that the data allow us to reliably estimate only the quarterly patient numbers of a doctor; we do not observe the maximal number on a given weekday, which would serve as a better proxy for the capacity. Capacity might depend on several nonobservable characteristics in the data, such as working hours, age, sex, number of assistants, infrastructure, or characteristics of the patient set (27). Capacity might also show seasonal variations and be lower during holiday periods. However, while these factors may certainly be helpful in better determining the capacity of a single doctor, they can be expected to have little impact on the overall systemic properties of the networks, such as the existence of a critical fraction of removed doctors on which nonresilient behavior sets in. The factors mentioned could certainly shift the position of but would not necessarily impact the comparison of individual districts.

There is a tradeoff between increasing resilience and decreasing overcapacity. One should be careful about interpreting the proposed resilience scores literally as excess capacity. Currently, we overestimate regional resilience by providing upper bounds for how much capacity is needed to avoid disruptions. For instance, if we choose such that removals of this size will lead to patients being lost with certainty, there is a chance that patients may already be lost for smaller shocks. We also assumed that patients have near-unlimited patience in searching for new doctors, when in reality, they may just abstain from consulting any further doctors after 1 or 2 unsuccessful attempts. Factors like these may decrease resilience and thereby, overcapacity. A diagnosis of overcapacity would also require an overrepresentation of providers with low financial income.

Our approach can be easily modified to include scenarios other than removals of doctors. There could be surges in patient numbers due to an influx of refugees (28) or an epidemic (29). The method can be transferred to other settings as long as the construction of a patient-sharing network is feasible (i.e., there is a negligible number of isolated nodes or groups of nodes). The structure of patient-sharing networks has already been studied in the United States (15–19), Canada (30), Italy (31), and Australia (32). To transfer our model to other settings, one would, therefore, need to 1) identify suitable data, 2) identify the relevant health sector (e.g., primary care), and 3) confirm that the networks are connected (no substantial isolated components). Most model parameters are estimated from historical data and therefore, take the heterogeneity of providers and health care delivery models explicitly into account. For instance, in Austria, each federal state has its own social security institution (as do certain occupation groups), which could confound the results. In the regression analysis, we showed that our regional resilience indicators are not driven by such state-level effects, while adjustments might be necessary to compare doctor-level results across different federal states. Finally, it should be noted that the underlying dataset is more than 10 y old and therefore, cannot be expected to adequately represent the current situation in Austria.

Our results clearly show that the resilience of health care systems cannot be described by trivial summary statistics, such as physician density. Resilience must be understood and measured as the property of how networked systems absorb and restructure themselves in response to shocks (5). We show how resilience can be quantified and used to aid decisions on optimal allocations and how investments for the increase of regional PCP densities would be most beneficial. We can estimate the systemic relevance of individual providers and therefore, identify which providers it would be particularly important to replace immediately on their retirement.