Stress-testing the Resilience of the Austrian Healthcare System Using Agent-Based Simulation

Patients do not access physicians at random but rather via naturally emerging networks of patient flows between them. As retirements, mass quarantines and absence due to sickness during pandemics, or other shocks thin out these networks, the system might be pushed closer to a tipping point where it loses its ability to deliver care to the population. Here we propose a data-driven framework to quantify the regional resilience to such shocks of primary and secondary care in Austria via an agent-based model. For each region and medical specialty we construct detailed patient-sharing networks from administrative data and stress-test these networks by removing increasing numbers of physicians from the system. This allows us to measure regional resilience indicators describing how many physicians can be removed from a certain area before individual patients won't be treated anymore. We find that such tipping points do indeed exist and that regions and medical specialties differ substantially in their resilience. These systemic differences can be related to indicators for individual physicians by quantifying how much their hypothetical removal would stress the system (risk score) or how much of the stress from the removal of other physicians they would be able to absorb (benefit score). Our stress-testing framework could enable health authorities to rapidly identify bottlenecks in access to care as well as to inspect these naturally emerging physician networks and how potential absences would impact them.


Introduction
Practising physicians form the backbone of our healthcare systems to provide the population with access to health care, i.e. "the timely use of personal health services to achieve the best health outcomes" 1 . Access to health care might be hampered by structural barriers including the number, type, concentration and location of care providers, but also by how quickly they can be accessed. 2 Consequently indicators that seek to quantify the access to health care are often derived from these metrics, such as the number of care providers per capita in a given region 3

. It was recently
shown that indicators such as the density of care providers fall short of capturing structural barriers to health care that arise from the fact that patients do not access physicians at random 4 . Rather, physicians are embedded in informal and emergent networks of flows of patients between them [5][6][7][8][9][10] .
These networks might, for instance, emerge because of geographic proximity, of one physician tending to recommend another one, or due to physicians acting as holiday substitutes to each other, and thereby encode how likely a patient will access a given care provider given that s/he has already accessed a specific other one. Such network-derived indicators for access to health care can paint completely different pictures on systemically important physicians with respect to quantitative indicators that solely focus on the density of care providers 4 .
The functionality of the healthcare system in many developed countries is increasingly challenged by demographic shifts both in the patient population 11 and the health workforce 12 (retirement waves), climate-related risks such as extreme weather 13, 14 as well as future and currently emerging infectious diseases 15 . Yet, surprisingly little attention is given to the problem of quantifying and stress-testing the resilience of national healthcare systems, as opposed to, e.g., economic or financial systems where such stress tests are a common risk management practice 16 . In this work we aim to fill this gap by developing a stress-testing framework to quantify the resilience of primary and secondary care by physicians in Austria using agent-based simulations. We adopt the notion of resilience as the ability to manage and regulate the adaptive capacities of a complex adaptive system that can be described via networks 17 . In particular, we consider how regional access to primary and secondary care by physicians depends on the capacities of individual physicians when confronted with shocks that reduce the number of available physicians (be it through retirement, quarantine, or other external shocks). We model the system's response to such a shock in terms of a restructuring of the emergent patient-sharing network and quantify how likely patients are to find new care providers with sufficient capacity within a specific location and time period. The agent-based model is fully calibrated to observational health care access data that covers approx. 100 million visits of 7,630,498 patients at 9,580 physicians of 13 different specialties in Austria.
Our approach can be briefly outlined as follows. Initially, all physicians are available. We consecutively remove individual physicians of a certain specialty and simulate the ensuing redistribution of patients within the healthcare system as patients of the removed physicians try to find a new health care provider. Patients choose their new physicians with probabilities that are proportional to the number of shared patients between the removed and the potential new physician. Based on their remaining free capacity, the newly chosen physicians will either accept some additional patients or reject them. The rejected patients will continue to look for a new available physician in the next time step, until they reach the maximum number of physicians they 4 are willing to contact (or a maximum distance they are willing to travel). If patients have not been accepted by a given number of physicians, these patients give up on getting an appointment and become 'lost'. We assume the patients' starting locations to be the municipality of their initial physician. The one-by-one removal process of physicians continues until all physicians become unavailable.
In this work we derive two sets of indicators to measure the regional resilience of primary and secondary care sectors. One set of indicators quantifies properties related to the resilience of the entire sector (e.g. specialty) in a given region, for example internal medicine in Vienna. The second set of indicators quantifies properties of individual physicians. For each physician we, therefore, define a risk and a benefit score. The risk score of a physician quantifies how much the removal of that physician affects other physician's capacities (the higher this score, the more likely it is that other physicians will exceed their capacity should the given physician be removed). The benefit score, on the other hand, relates to the free capacity of a physician.
By considering the percentage of physicians with a given specialty in a region that have a risk or benefit score higher than the nation-wide average, we compute regional risk and benefit levels for each specialty. Furthermore, for each region and specialty we measure the free capacity (defined as the percentage of physicians that can be removed before 20% of the initially free capacity has been filled up) and an indicator for the number of lost patients (percentage of physicians to be removed before 1% of the patients in a region don't find a new care provider in the model).
Agent-based models often produce complex and high-dimensional outputs that can be hard to 5 interpret for both technical and non-technical experts 18 . To address this challenge we develop a visualisation strategy for the indicators; see Figure 1. There we show the state-and specialty-specific resilience indicators that are described by glyphs of different shapes and colour. State-specific results include aggregated values of four indicators and are depicted by a diamond shaped glyph within a circle (a). The height and width of the diamond shape represent the states' resilience in terms of free capacity and lost patients, respectively. The colour-filled arcs within the semicircles describe the state level percentages of physicians that have higher than state-average risk and benefit scores. The physician-or specialty-specific resilience is described by two of the previous indicators, namely risk and benefit scores (c). These are visualised by the colour-coded halves of a heart-shaped symbol that can represent results of individual physicians as well as aggregated results of a medical specialty.
To illustrate how we measure these indicators in the network model, see Fig. 1 (c,d). There we show the network neighbourhood of a physician ("physician 1") that becomes unavailable and has a relatively high risk score, meaning that his/her removal will affect many other physicians due to an inflow of new patients (c). As a result, physician 2 with a relatively low benefit score (meaning low capacity to accept new patients) becomes unavailable, and some patients need to contact yet another physician to find a new care provider and may potentially get lost (d).
6 Figure 1: Schematic overview. a) Custom glyph used to describe state-level results. b) Custom glyph for aggregated specialty-or physician-specific resilience indicators. c) and d) Example of physicians' patient sharing network and patient displacement: (c) physician 1 becomes unavailable and a share of their patients try to re-locate to physician 2. As a result, physician 2 becomes unavailable as well (d). (e) State-specific result panel for selected medical specialties (psychiatrists) including averaged resilience indicators, numbers of patients and physicians and distributions of risk and benefit scores.

Agent based model
Data-sets To construct the patient sharing network, we combine two data-sets: patient contacts for Austrian physicians as derived from administrative data, and a data set of physician's opening hours that was scraped from a health care information platform (www.herold.at) in March 2020 1 . Some physicians practice more than one medical specialty and/or practice in more than one municipality. Therefore, in the following we always refer to a unique (physician, specialty, municipality) combination, when talking about the number of physicians. A physician who is both, say, a general practitioner and an internist, or a physician who operates as a general practitioner in two municipalities is counted as two different physicians.

Patient contacts
The health record data provided by the Austrian Federal Ministry of Health contains all recorded patient visits (approx. 103,000,000) over a 1-year period from January 1, 2018 to December 31, 2018. The data contains patients' contacts to physicians from 13 different medical specialties and includes a pseudonymised unique ID for every patient and physician, the date of contact as well as the locations of physicians in form of 5-digit municipality IDs. The first digit of the municipality ID is further used to determine the federal state (9 states in total). The data set contains 9,580 unique combinations of (physician, specialisation, municipality). 1 Permission for scraping the data was obtained from the owners of the platform Patient sharing network To model connections between physicians, we use the patient contact data to construct a patient sharing network. In this network, every node is a physician. The network is described by the adjacency matrix A i,j , where entries a i,j describe the number of shared patients between physicians i and j, i.e. the number of patients that visited both physician i and physician j in a given period of time. The adjacency matrix is therefore symmetrical by definition and the network described by the adjacency matrix is undirected. To construct the adjacency matrix, we start with a patient's visit to physician i and consider all visits of the given patient to other physicians during a period of three months before and after the initial visit to physician i. All thus identified physicians are therefore sharing one patient with i. This process is repeated for all patient visits of physician i. If there is at least one shared patient between i and a given other physician j, there is an edge a i,j between i and j. The edge is weighted by the sum of shared patients. This procedure is repeated for all physicians, resulting in a weighted undirected patient sharing network that is represented by the adjacency matrix A i,j . To reduce spurious relations, the adjacency matrix can be thresholded by requiring a minimum number of shared patients and a maximum geographical distance between physicians i and j. The distance is measured as the direct distance between i's and j's municipalities. The location of municipalities was determined as either the geolocation of the affiliated municipal office or as the geolocation of the area with highest population density in the municipality using a population density map. If the distance between physicians i and j exceeds threshold d, the edge is removed and a i,j = a j,i = 0. Similarly, when the number of shared patients between i and j is smaller than threshold p, the edge is removed.
To model the patient sharing network within a medical specialty, only the subset of nodes 9 corresponding to physicians of the given specialty and their corresponding edges are used. Edges between physicians of different specialties are removed.
Capacity estimation To estimate a physician's maximum capacity, C, we aggregate the patient contact data with information about opening hours; see SI Section "Matching of patient contact data to opening hours" for details. For robustness, we provide two different estimates of capacity: one based purely on the number of patient contacts (C p ) and one based on the opening hours and patient contacts (C h ).
For the patient-contact-based capacity estimate C p , we calculate the median patient capacity (number of quarterly patient contacts) of the top 10% of physicians for each specialty. We then assign this median capacity to the rest of the physicians (not top 10% capacity), regardless of opening hours. Since this approach does not differentiate between physicians with different opening hours, it is prone to overestimating their capacity.
For the opening-hour-based capacity estimate C h , we divide physicians of every specialty into bins based on their total number of weekly opening hours with a bin size of 5 opening hours. For each bin, we then identify the top 10% of physicians in terms of their capacity (number of quarterly patient contacts) and calculate their median capacity. This median capacity is then assigned to all remaining physicians in the bin as their maximum patient capacity.
By merging the information of the two data-sets we get a comprehensive description of the majority of physicians in the country. We can initially describe each physician with specific characteristics: the specialty, municipality, number of patients seen per quarter (N ) and total capacity based on opening-hours (C h ). While the specialty and location are constant, the current number of patients and the remaining free capacity change during simulations. General descriptive information is shown in Tab. 1. The vast majority of physicians are general practitioners (n = 4, 967), while psychiatrists (n = 161) form the smallest specialty group. Patient variables in the model include their starting location (municipality) based on the physician they start at, the number of displacement steps they have already undergone, step-wise travelled distance and total travelled distance.

Model parameters and robustness
The model contains a number of parameters that can be varied and tested for robustness. We define a baseline parameter setting as follows: The maximum number of displacement steps before a patient is considered lost (s = 10), the median patient capacity of top c% of most visited physicians of certain specialty that is assigned to all other physicians as their maximum capacity (c = 10), the maximum distance in kilometres patients are willing to travel to a new physician (d = 100), the number of times the simulation attempts to find a physician within this distance from the starting location of a patient before choosing a physician located further away (i = 10), the minimum number of shared patients in the adjacency matrix for a valid connection   Maximum capacities C h were calculated based on opening hour information, using the top 10% of physicians per specialty.

Free capacity & lost patients
The initial free capacity in each federal state is reduced in each simulation step as more and more physicians become unavailable. We track how the remaining relative free capacity in each state is filled up/diminished until there is no capacity left and all the physicians are removed from the system. We define a critical limit of 20% of the remaining free capacity per state and track the average relative number of physicians that needs to be removed to reach this critical free capacity limit, L F C .
Similarly, we sum the number of lost patients and calculate the average relative number of cumulative lost patients in each federal state. The more physicians are unavailable, the lower the remaining free capacity and the higher the probability for a patient to get rejected at a new physician -after all physicians are removed, 100% of patients are lost. Analogous to the free capacity we define a critical lost patient limit of 1% and track the relative number of physicians removed before this point is reached, L LP .
Risk & benefit scores Based on the initial physicians' characteristics we can describe the negative and positive influence that physician contributes to the healthcare system by defining a risk and benefit score for each physician. The risk score R i represents physician i's risk and is measured by the average extra load of patients that their first-degree neighbours in the patient sharing network must bear in case of i's unavailability. We therefore define the risk score as where N i and N j stand for the number of patients of physician i and their neighbours j, C j is j's total capacity and w j describes the normalised connection weight between j and i. The range of the 13 risk score is [0, 1], where 0 corresponds to the lowest and 1 to the highest risk. The benefit score B i represents physician i's beneficial contribution to the healthcare system in terms of their initial free capacity. The benefit sores are normalised over all physicians of a given medical specialty and range as well [0, 1], with 0 being the lowest benefit and 1 the highest benefit.

Results
Regional resilience indicators To investigate indicators of regional resilience of the healthcare system, we report the development of the relative lost patients as well as the remaining free capacity as physicians are gradually removed in the simulation. We report averaged results over 100 simulation iterations. Figure 2   The slopes of this linear dependence, however, vary widely across federal states and relate to the lost patients indicator: the faster capacity decreases, the sooner patients get lost. As the system gradually approaches a state with no remaining capacity, there is no more free capacity to accommodate displaced patients and an increasing number of patients get lost.
In the SI, Figures S4 to S15, we show lost patients and free capacity as a function of the percentage of removed physicians for all other specialties. A summary for all states and specialties and the relative number of physicians that can be removed from the system until critical limits are reached (20% for the remaining free capacity and 1% for lost patients) is shown in Fig. 3.
The resilience indicators vary widely across federal states and specialties. Resilience indicators in the state of Styria are typically on the lower side, whereas structures in Vorarlberg appear to be substantially more resilient. Psychiatrists, neurologists, surgery, internal medicine, orthopaedics and radiology show a tendency to more resilient care networks compared to ophtalmology, dermatology, or urology. This can also be seen in Fig. 4, where we show nationwide averages of these resilience indicators.
In addition to regional resilience of the healthcare system as a whole, we also investigate the individual risk and benefit a given physician contributes to the system. To this end, for each specialty j we average the individual risk(risk j ) and benefit (benefit j ) scores over all federal states to obtain nationwide results. Figure 4 (b) shows the mean and standard deviation for each specialty averaged over all federal states.
In general, there is a tendency that specialties with a high risk score and/or a low benefit score are also less resilient in terms of lost patients and free capacity. To quantify this observation, we evaluate a linear regression model of the form L F C j (L LP j ) ∼ r j · risk j + b j · benefit j + const., where L F C j and L LP j are the state-averaged percentages of physicians of specialty j removed before reaching critical limits of free capacity and lost patients, respectively. The averaged effect sizes and standard deviations of coefficients r and b over all j confirm the correlation between these indicators as r = −11(±10), b = 25(±28) for L F C and r = −70(±35), b = 58(±86) for L LP . However, there is no discernible correlation between the risk and benefit scores across the specialties. This suggests that these scores capture two structurally independent properties that together determine the resilience of a specialty. This also motivates the use of the individual-level risk and benefit scores to assess the contributions of physicians to the systemic resilience of these regional care sectors. Darker colours indicate higher number of physicians that can be removed before exceeding critical limits.
Tables including the standard deviation of simulation results are given in Tables S2 and S3. Robustness We performed several robustness tests to assess the sensitivity of our results. We find that varying the maximum capacity parameter c (c = 20, 30, 40) introduced most variation in the results. However, this variation does not change the regional resilience indicators qualitatively, meaning that the specific values for overstepping critical limits of lost patients or free capacity will differ for states when using different c, but it will not convert a "more resilient" state to a "less resilient" state compared to others (ranking remains robust). For the purpose of this stress-test we selected c = 10 as it necessitates a large proportion of physicians to increase their available capacities.
Increasing the probability to approach physicians at random instead of following the links in the patient-sharing network, i.e. setting α = 0.15, showed the second largest variation in simulation outcomes. Again, rankings remained robust as increasing α gradually washes out the network effects. In the absence of such network effects (for α = 1) our resilience indicators would be determined by physician density alone. Other simulation parameters only show minor changes in simulation results. Every other model parameter than c and α was tested with at least one additional setting (s = 5, d = 50, i = 5, p = 1) to assess robustness of the results.

Discussion
Fast progress in digital technology stimulates debates on new consultation methods. Especially the Covid-19 pandemic encouraged and increased the implementation of online consultations in many fields of medicine 19,20 . However, patients in regions with more deprived populations appear to be less likely to benefit from telehealth 21 . Patients with inadequate digital devices or in regions with poor network connectivity will have problems to access telehealth, possibly even widening already existing treatment gaps. On the other hand, remote GP but also specialist consultations may be helpful and convenient for many patients with chronic diseases in need of regular follow-up visits.
However, a certain number of face-to-face consultations will remain necessary to provide adequate patients care in the future. Therefore, although telehealth may help to diminish problems of physicians' shortage in some areas in the future, allowing continued access to services, we will still need a sufficient number of experienced physicians with capacities to perform online consultations in addition to standard clinical care.
Experiences during the Omicron-variant-related surge in the SARS-CoV-2 pandemic have demonstrated that scenarios in which ten to twenty percent of the health workforce are absent over several weeks can indeed happen 22 . As health care demand is further expected to increase in developed countries as the population is ageing and the prevalence of chronic conditions is increasing 11 , there is an increasing need to know how far the available physician supply can be stretched to still deliver the desired quality of care. In this work we proposed novel regional resilience indicators to quantify the point beyond which further decreases in physician densities would begin to severely impair the population's ability to access primary or secondary care. We find strong heterogeneities in these resilience indicators across different medical fields and regions.
Our findings show that the naturally emerging patient sharing networks that can be reconstructed from routinely collected administrative data are a valuable source of information when assessing the resilience of the healthcare system.

Our resilience indicators can be disaggregated into indicators for individual physicians to
quantify the extent to which their hypothetical removal would stress the system (risk score) or how much of the stress from the removal of other physicians they would be able to absorb (benefit score).
These indicators allow health authorities to rapidly identify physicians with sufficient capacities to take over supply and alert them that they might receive larger than expected inflows of patients in the situation where other physicians retire. This could help planners to proactively anticipate whether such retirements could cause bottlenecks in certain regions and whether the now vacant position should be replaced with a higher level of priority. To highlight this potential of our approach, we have developed an interactive visualisation framework that allows health care planners and professionals to browse and intuitively access the regional and individual resilience indicators as well as to inspect these naturally emerging physician networks and how potential absences would impact them.
Our study has several limitations. Firstly, we focused on patient-sharing networks among physicians of only the same medical specialty. However, especially in rural areas GPs are used to cover and treat a broad spectrum of diseases and they may be able to attenuate the gap caused by shortage of some specialists 23 . On the other hand, specialists in internal medicine could temporarily replace GPs in some areas. Moreover, in Austria many tasks of primary care are still transferred or covered by outpatient departments of bigger hospitals in the cities. Therefore, shortfall of physicians might have different effects in rural and urban areas.
Secondly, the results of the agent-based simulation depend on the assumptions regarding various model parameters. We performed several robustness tests to better understand the sensitivity of our results. We find that our results are most sensitive to changes in the maximum capacity 21 parameter c and teleportation parameter α. However, we observed no substantial changes in the rankings of states and specialties, respectively.
Finally, the sex of physicians should be included in future analyses because recent trends show that the number of female physicians is steadily increasing but they appear to have different demands and preferences regarding work-life balance, working models and fields of activities and duties 24 . Moreover, they need specific infrastructures like family-friendly workplace.
Our work demonstrates that the resilience of primary and specialist care is an emergent property of formal and informal networks that physicians form amongst themselves. While fiscal constraints and concerns regarding an ageing physician workforce in rural areas are growing worldwide 26 , we find that these care networks indeed possess tipping points in terms of their ability to provide care to the entire population. Removing or losing an alleged excess capacity or oversupply of physicians might inadvertently push healthcare systems closer to their tipping points.
Using our indicators allows health authorities to quantify how close they are to these tipping points and thereby strike a more data-informed balance between system resilience and effectiveness. on a given day are assumed to be open for 2 hours on that day. If there is more than one entry of the same physician but with different opening hours (29 entries in total), we keep the entry with the longer opening hours. We then match the 29 specialities given in the opening hours data set (see table S1) to the 13 specialities given in the patient contact data (see Table 1). We drop 9 specialities -such as medical and chemical laboratory diagnostics -that are not corresponding to any of the specialities in the patient contact data (90 physicians). This leaves us with a total of 6,604 unique physicians.
We then probabilistically assign the opening hour data to the patient contact data follow-27 ing a two-step process: We first look for all instances where there is only a single physician with a given speciality in a given district in both data sets and match these physicians. This We then assign physicians with a given speciality and within a given municipality based on the difference between their opening hours and patient contacts by minimizing We also introduce a threshold = 0. For the unmatched physicians, we impute the opening hour information in the following way: for every unmatched physician, we look for all successfully matched physicians with the same specialisation and similar number of patients. We first look for all physicians that have the same number of annual patients ±10. If no physicians are found, we increase the search interval to ±100 or ±1000 annual patients. We then calculate the average opening hours for each day of the week 29 for these successfully matched physicians and use this information to impute the missing opening hour information for the unmatched physicians.

B Progress of resilience indicators in states
The following figures (from Figure S4 to Figure Fig. 3 in the main text.

D Interactive online visualisation tool
The tool contains three components (see Fig. S16 and S17): 1) an overview of resilience indicators per federal state and medical specialty; 2) a detailed view of indicators for the selected specialty in the selected federal state, including risk & benefit scores, information on free capacity & lost patients; 3) the physicians' network within a medical field, displaying indicators per physician. In order to allow the comparison of resilience indicator on an aggregate level, we developed two custom glyphs that respectively summarise the characteristics of each federal state and medical specialty. To encode the aggregate resilience of a medical specialty, we present the average values of the risk and benefit scores of physicians in a field within the two halves of a heart shaped glyph depicted in Fig. S19. We use a similar approach as in the semicircle representation for federal states: the colour of the left half of heart shows the risk score and the right half shows the benefit score. The overview panel of resilience indicators also serves as a filter for selecting a particular 46 medical field in a particular federal state. After the selection, three sub-components are provided: First, a basic profile that contains four attributes on state level (free capacity, lost patients, risk and benefit levels) and two attributes on specialty level (risk and benefit scores) that are described as an aggregate as well as by their distributions, and finally the number of physicians and patients. The second component shows the model results during a physician removal step (Fig. S20). The third component displays the list of physicians in the selected specialty.
The model results contain three charts: remaining free capacity, average displacement steps and average travel distance. The charts visualise the value changes in respect to the percentage of removed physicians. Users can use a slider to adjust the value from "0%" to "100%", in order to remove a certain amount of physicians (on national level) and compare the values. The physicians' network component is used to illustrate the displacement steps of a "lost" patient. It assumes that the current chosen physician always turns out to be unavailable, so the user has to choose another physician based on the current physician's network (see Fig. S22). We expect the user can learn the fundamentals of this simulation system better by interacting with a visual network.
In the visit history in Fig. S23 on the top of the panel, users can view the steps they have made and undo the previous step. Since in our model we also assume that patients can (and are willing to) travel only a limited distance from their starting location to find a new physician and is 48 Figure  only willing to contact a limited number of new physicians, the user cannot select a physician who is more than 100 km away from the starting location of their initial physician, or over 10 steps away in the physician network. As such, the user always plays the role of a "lost" patient in order to learn how patients can get lost in the system.