The first implication is extending the multisystemic perspectives in resilience development by seeing resilience as a complex adjustment process involving multiple systems and contexts. The risk of crises accumulates as the number of disturbance events rises beyond the resilience level. Abrupt, significant, temporary, and frequent changes in surrounding organizations may cause unintended difficulties for individuals, primarily when changes concern essential resources or support for one's immediate stabilization and grounding. For people who struggle with acute abruptions and disturbances, these changes may become another burdensome object that people need to become acquainted with quickly, which is likely to contribute to emotional disturbance in the sense of experiencing the unknown, uncertainty, and feelings of powerlessness, helplessness, and The first implication is extending the multisystemic perspectives in resilience development by seeing resilience as a complex adjustment process involving multiple systems and contexts. The risk of crises accumulates as the number of disturbance events rises beyond the resilience level. Abrupt, significant, temporary, and frequent changes in surrounding organizations may cause unintended difficulties for individuals, primarily when changes concern essential resources or support for one's immediate stabilization and grounding. For people who struggle with acute abruptions and disturbances, these changes may become another burdensome object that people need to become acquainted with quickly, which is likely to contribute to emotional disturbance in the sense of experiencing the unknown, uncertainty, and feelings of powerlessness, helplessness, and hopelessness.

hopelessness. The second implication relates to strengthening and reinforcing protective factors in producing systemic efforts of strategy development and policy design in population resilience. At an individual level, resilience developed through recovering from previous shocks of the exact nature is critical as it provides higher tolerance of distress and adversity. The recovery process requires individuals to constantly build on personal assets such as self-esteem and emotional positivity. It is also critical to develop adequate insights into one's external reality, nuanced emotions, and intrapsychic experiences. Individuals can The second implication relates to strengthening and reinforcing protective factors in producing systemic efforts of strategy development and policy design in population resilience. At an individual level, resilience developed through recovering from previous shocks of the exact nature is critical as it provides higher tolerance of distress and adversity. The recovery process requires individuals to constantly build on personal assets such as self-esteem and emotional positivity. It is also critical to develop adequate insights into one's external reality, nuanced emotions, and intrapsychic experiences. Individuals can present with psychological flexibility, reflect, and wonder about goals and visions of oneself and life, as well as present with adaptive coping skills and strategies, which support

individuals' counter against appraisal of negative emotions and destructive responses. Involving multi-sector organizations such as the education, workplace, and public service sectors across systems is essential, as resilience development is a process of experiencing human relationships, connections, rapport, and trust across the boundaries of multiple systems. Resilience of multiple systems serves not only certain designed functions, but also forms stable systems to protect individuals, families, and even groups from being shattered by crises and forceful interruptions. As the model suggested, systemic efforts such as equipping parents with skills and knowledge about crises and providing supportive resources, strategic guidelines, and financial and employment security are needed to facilitate the growth of resilience.

The third implication is conceptualizing organizational resilience from individual resilience development perspectives. Persistent crises require organizations to respond with mitigations, but short-term adjustments that only focus on securing organizational performance and functioning without considering individual resilience can increase the accumulative disturbances individuals face within the short term, which hinders the development of psychological resilience during crises. Organizational and individual learnings are critical to help the growth of resilience and decrease the number of people becoming increasingly stressed during a crisis. When crises and shocks to the system regard health and well-being, such as COVID-19, measurements of organizational resilience should expand from organizational functionality and performance. An individual–environment view of organizational resilience should incorporate indicators such as to what extent the organization provides insights and learnings from the adverse events for individuals, and to what extent the organization's capacity to adapt to disturbances and mobilize resources to sustain changes and provide transformational adjustments.

The last implication regards COVID-19 lessons. Significant changes have been observed in conducting and receiving health and well-being services. As all parties strive to adjust to the long-pandemic, we wonder how we may learn and reflect from our experiences of the pandemic and continuously support people in need. When the population faces constant and enduring waves of a pandemic, policies and multiple sectors should facilitate the development of psychological assets such as self-esteem, psychological flexibility, adaptive self-care or coping strategies, and supportive social welfare. The present reality is that the pandemic has a high level of uncertainty, and the unknown may not be eliminated, which requires us to live with these disturbing and uncomfortable dynamics. A resilience perspective can cultivate tolerance of disturbances and distress, reinforcing the population's growth to strive for a happy life.

In terms of limitations, while we used the COVID-19 pandemic as an example to explore the interconnections between individual and organizational resilience, the model is not calibrated with empirical data and does not consider variations across subpopulations. Emerging longitudinal research shows that young generations aged 19–30 years and females had a higher risk of distress during the pandemic and its lockdowns [57]. Additionally, younger age (<40 years), female gender, psychological illness, student status, exposure to social media/news, and unemployment are common risk factors that have been shown to be associated with mental distress caused by the pandemic [18]. The model was not calibrated with empirical data, as the focus of this model is to provide theoretical exploration. Empirical evidence decides the realism and reliability of the model's validity, which requires both structural and behavior-over-time data [58]. While we believe theory-based modeling was useful and proper for the purpose of this paper, the model was not calibrated with behavior-over-time data; thus, the model should be considered exploratory with tentative and uncertain conclusions. To improve the evidence level of the model, future studies can calibrate the model using data on the population's health and well-being, service provision, and environmental changes. Future research can further include the impact of accumulative risk on different socioeconomic or age groups to show the different levels of vulnerability.

The model provided theoretical insights but was overall simplified at quite a high level. A more systematic review of theoretical orientations or models in clinical psychology can be used to improve the structural conceptualization. For example, resilience can be closely linked to psychological development, insights, judgments, resourcefulness, behavioral regulation, distress tolerance, etc. For organizational resilience, we simplified the conceptualization of organization resilience by using the number of staff in the health and well-being service sector. In practice, different agencies have broader approaches to health services, such as crisis intervention and community-based health services. Agencies that provide community health services and care might have more significant challenges in responding to crisis interventions during multiple pandemic waves. Future research can integrate broader service provision challenges by different agencies. Meanwhile, it could be meaningful to include resilience in different cultures, communities, and groups, so this concept and thinking can be further nuanced to different cultural norms and identities.

**Author Contributions:** Conceptualization, K.Z.; methodology, K.Z.; modeling, K.Z.; writing–-original draft preparation, K.Z. and M.Z.; writing—review and editing, K.Z. and M.Z.; visualization, K.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** K.Z. was funded by the National Institute for Health and Care Research (NIHR) School for Public Health Research (SPHR), grant reference number PD-SPH-2015. The views expressed are those of the authors and not necessarily those of the NIHR or the Department of Health and Social Care.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors appreciate Christina Gkini, Jefferson Rajah, and Birgit Kopainsky for their valuable feedback on the model. They are also grateful for feedback received at the 2022 System Dynamics Conference in Frankfurt, Germany.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A Model Equations**

\*\*\*\*\*\*\*\*\*\* Crises\_multiple\_shocks: \*\*\*\*\*\*\*\*\*\* First\_shock\_duration = 5 UNITS: month First\_shock\_intensity = 3 UNITS: dmnl First\_shock\_start = 20 UNITS: month First\_shock\_stop = First\_shock\_start+First\_shock\_duration UNITS: month Interval\_to\_second\_shock\_start = 6 UNITS: month Interval\_to\_third\_shock\_start = 6 UNITS: month Multiple\_Shocks[First\_wave] = STEP(Risk\_intensity\_of\_each\_shock[First\_wave], First\_shock\_start)-STEP(Risk\_intensity\_of\_each\_shock[First\_wave], First\_shock\_stop) UNITS: dmnl Multiple\_Shocks[Second\_wave] = STEP(Risk\_intensity\_of\_each\_shock[Second\_wave], Second\_shock\_Start)-STEP(Risk\_intensity\_of\_each\_shock[Second\_wave], Second\_shock\_stop) UNITS: dmnl Multiple\_Shocks[Third\_wave] = STEP(Risk\_intensity\_of\_each\_shock[Third\_wave], Third\_shock\_start)-STEP(Risk\_intensity\_of\_each\_shock[Third\_wave], Third\_shock\_stop)

UNITS: dmnl "MULTIPLE\_SHOCKS\_Switch\_1=\_crisis\_on" = 1 UNITS: dmnl Risk\_intensity\_of\_each\_shock[First\_wave] = First\_shock\_intensity UNITS: dmnl Risk\_intensity\_of\_each\_shock[Second\_wave] = Second\_shock\_intensity UNITS: dmnl Risk\_intensity\_of\_each\_shock[Third\_wave] = Third\_shock\_intensity UNITS: dmnl Second\_shock\_duration = 5 UNITS: month Second\_shock\_intensity = 3 UNITS: dmnl Second\_shock\_Start = Interval\_to\_second\_shock\_start+First\_shock\_stop UNITS: month Second\_shock\_stop = Second\_shock\_duration+Second\_shock\_Start UNITS: month Third\_shock\_duration = 5 UNITS: month Third\_shock\_intensity = 3 UNITS: dmnl Third\_shock\_start = Second\_shock\_stop+Interval\_to\_third\_shock\_start UNITS: month Third\_shock\_stop = Third\_shock\_start+Third\_shock\_duration UNITS: month \*\*\*\*\*\*\*\*\*\* Disturbance\_and\_Psychological\_resilience: \*\*\*\*\*\*\*\*\*\* Disturbance\_events\_that\_can\_be\_dealt\_with\_resiliency = Psychological\_Resilience// Resilience\_needed\_per\_disturbance\_event UNITS: event Fractional\_emotional\_tolerance\_acquired\_from\_recovering = IF P3\_Collective\_growth =1 THEN Normal\_fractional\_emotional\_tolerance\_acquired\_from\_recovering\*Intensity\_of \_P3\_on\_distress\_tolerance ELSE Normal\_fractional\_emotional\_tolerance\_acquired\_from \_recovering UNITS: dmnl Goal\_of\_resilience = Resilience\_needed\_per\_disturbance\_event\*Disturbance\_events UNITS: resilience Initial\_individual\_resilience = INIT(Goal\_of\_resilience) UNITS: resilience Normal\_fractional\_emotional\_tolerance\_acquired\_from\_recovering = 1 UNITS: dmnl Psychological\_Resilience(t) = Psychological\_Resilience(t - dt) + (Resilience\_adjusting \_by\_Recovering\_from\_shocks) ∗ dt INIT Psychological\_Resilience = Initial\_individual\_resilience UNITS: resilience Psychological\_resilience\_ratio = EXP(Psychological\_Resilience)/(1+EXP(Psychological\_Resilience)) UNITS: dmnl Resilience\_adjusting\_by\_Recovering\_from\_shocks = (Goal\_of\_resilience-Psychological\_Resilience)\*(Fractional\_effect\_of\_tolerating\_emotional\_distress\_through\_recovering \_process\_per\_month)\*Environment\_stabilization\_ratio UNITS: resilience/month Resilience\_needed\_per\_disturbance\_event = 0.053

UNITS: resilience/event \*\*\*\*\*\*\*\*\*\* Environment\_adjustments: \*\*\*\*\*\*\*\*\*\* Adjustments\_Gap = MAX(0, SMTH3(Total\_adjustments\_needed-Environment\_Stabilized\_Adjustments, Organization\_response\_time)) UNITS: adjustment Environment\_adjusting = MAX(0, (Adjustments\_Gap)//Time\_needed\_to\_develop\_ temporary\_adjustments) UNITS: adjustment/month Environment\_Stabilized\_Adjustments(t) = Environment\_Stabilized\_Adjustments(t dt) + (Environment\_stabilizing - "Re-adjusting") ∗ dt INIT Environment\_Stabilized\_Adjustments = Initial\_stabilised\_adjustments UNITS: adjustment Environment\_stabilizing = Fractional\_adjustments\_that\_move\_to\_stabilized\_adjustments\* Environment\_Temporary\_Adjustments//Time\_needed\_to\_settle\_permanent\_solutions UNITS: adjustment/month Environment\_Temporary\_Adjustments(t) = Environment\_Temporary\_Adjustments(t dt) + (Environment\_adjusting + "Re-adjusting" - Environment\_stabilizing) ∗ dt INIT Environment\_Temporary\_Adjustments = Initial\_temporary\_adjustments UNITS: adjustment Fractional\_adjustments\_that\_move\_to\_stabilized\_adjustments = 0.3 UNITS: dmnl Initial\_stabilised\_adjustments = Total\_adjustments\_needed UNITS: adjustment Initial\_temporary\_adjustments = Total\_adjustments\_needed\*Time\_needed\_to\_settle\_ permanent\_solutions/(Time\_needed\_for\_revisiting\_adjustments\*Fractional\_adjustments\_ that\_move\_to\_stabilized\_adjustments) UNITS: adjustment Normal\_organizational\_response\_time = 3 UNITS: month Normal\_time\_needed\_to\_develop\_temporary\_adjustments = 4 UNITS: month Normal\_time\_needed\_to\_revisit\_adjustments = 12 UNITS: month Normal\_time\_needed\_to\_settle\_permanent\_adjustments = 12 UNITS: month Organization\_response\_time = IF P1:\_Environment\_Fast\_Adaptation=1 THEN Normal\_organizational\_response\_time/Intensity\_of\_P1\_on\_organizational\_response\_time ELSE Normal\_organizational\_response\_time UNITS: month "Re-adjusting" = Environment\_Stabilized\_Adjustments//Time\_needed\_for\_revisiting\_adjustments UNITS: adjustment/month Time\_needed\_for\_revisiting\_adjustments = IF P1:\_Environment\_Fast\_Adaptation=1 THEN Normal\_time\_needed\_to\_revisit\_adjustments/Intensity\_of\_P1\_on\_revisiting\_adjustments ELSE Normal\_time\_needed\_to\_revisit\_adjust ments UNITS: month Time\_needed\_to\_develop\_temporary\_adjustments = IF P1:\_Environment\_Fast\_Adaptation=1 THEN Normal\_time\_needed\_to\_develop\_temporary\_adjustments/Intensity\_of\_- P1\_on\_temporary\_adjustment\_time ELSE Normal\_time\_needed\_to\_develop\_temporary- \_adjustments

UNITS: month

Time\_needed\_to\_settle\_permanent\_solutions = IF P3\_Collective\_growth=1 THEN Normal\_time\_needed\_to\_settle\_permanent\_adjustments/Intensity\_of\_P3\_on\_stabilised- \_adjustment ELSE Normal\_time\_needed\_to\_settle\_permanent\_adjustments UNITS: month \*\*\*\*\*\*\*\*\*\* Health and Well-being\_service: \*\*\*\*\*\*\*\*\*\* Changes\_of\_employee = CAPACITY\_RESPOND\_STWITCH∗ (Indicated\_demanded- \_number\_of\_MH\_staff-HW\_service\_staff)//Time\_to\_hire UNITS: people/month HW\_service\_staff(t) = HW\_service\_staff(t - dt) + (changes\_of\_employee) ∗ dt INIT HW\_service\_staff = Initial\_MH\_service\_staff UNITS: people Indicated\_demanded\_number\_of\_MH\_staff = Total\_sessions\_demand/Number\_of- \_sessions\_per\_staff\_per\_month//Target\_Waiting\_Time UNITS: people Normal\_target\_waiting\_time = 4 UNITS: month Normal\_time\_to\_hire = 12 UNITS: month Number\_of\_sessions\_per\_staff\_per\_month = 100 UNITS: session/people/month Remaining\_monthly\_sessions\_capacity = MAX(0, Total\_monthly\_session\_capacity-Total\_monthly\_sessions\_occupied) UNITS: session/month Target\_Waiting\_Time = IF P2\_Mental\_Wellbeing\_Service\_Fast\_Response=1 THEN Normal\_target\_ waiting\_time/Intensity\_of\_P2\_on\_waiting\_time ELSE Normal\_target\_waiting\_time UNITS: month Time\_to\_hire = IF P2\_Mental\_Wellbeing\_Service\_Fast\_Response=1 THEN Normal \_time\_to\_hire/Intensity\_of\_P2\_on\_time\_to\_hire ELSE Normal\_time\_to\_hire UNITS: month Total\_monthly\_session\_capacity = HW\_service\_staff\*Number\_of\_sessions\_per\_staff \_per\_month UNITS: session/month Total\_monthly\_sessions\_occupied = "Highly stressed\_people\_with\_HW\_services" \*Frequency\_of\_sessions\_attended\_per\_month\_per\_people UNITS: session/month Total\_number\_of\_service\_sessions\_needed\_per\_person = 10 UNITS: session/people Total\_sessions\_demand = ("Highly\_stressed\_people\_without\_HW\_services"+"Highly \_stressed\_people\_with\_HW\_services")\*Total\_number\_of\_service\_sessions\_needed\_per \_person UNITS: session \*\*\*\*\*\*\*\*\*\* Initial\_numbers: \*\*\*\*\*\*\*\*\*\* Fractional\_recovering\_by\_own = Fraction\_of\_self\_recover/"Time\_to\_self-recover" UNITS: 1/month Fractional\_Stressing\_up = Normal\_fraction\_of\_highly\_stressed\_symptoms/Time\_to \_change\_stress\_level UNITS: 1/month Fractional\_using\_service = Fraction\_of\_service\_using\_among\_highly\_stressed\_ people/Target\_Waiting\_Time

UNITS: 1/month INIT\_highly\_stressed\_not\_using\_HW\_services = INIT\_low\_stress\_population\*(fractional\_Stressing\_up)/(fractional\_recovering\_by\_own+Fractional\_using\_service) UNITS: people INIT\_highly\_stressed\_using\_HW\_services = INIT\_highly\_stressed\_not\_using\_HW \_services\*Fractional\_using\_service\*Time\_to\_recover/(1+Time\_to\_recover\*Fraction\_of\_ dropping\_out\_per\_month) UNITS: people INIT\_low\_stress\_population = 10000 UNITS: people Initial\_MH\_service\_staff = INIT(Indicated\_demanded\_number\_of\_MH\_staff) UNITS: people Sum\_of\_population\_stocks = "Low/mild\_stress\_people" + "Highly-stressed\_people \_with\_HW\_services" + "Highly-stressed\_people\_without\_HW\_services" UNITS: people \*\*\*\*\*\*\*\*\*\* Multiple\_shocks\_and\_disturbance: \*\*\*\*\*\*\*\*\*\* Decreasing\_disturbance\_events = (Effect\_of\_environmental\_resilience\_on\_processing \_disturbance\_events\*Effect\_of\_individual\_resilience\_on\_processing\_disturbances)\*Disturbance\_events/Time\_needed\_to\_resolve\_events UNITS: event/month Disturbance\_events(t) = Disturbance\_events(t - dt) + (Increasing\_disturbance\_events - Decreasing\_disturbance\_events) ∗ dt INIT Disturbance\_events = Initial\_disturbance\_events UNITS: event Disturbance\_to\_temporary\_adjustments\_ratio = Environment\_Temporary\_Adjustments//Max\_adjustment\_depending\_on\_resiliency UNITS: dmnl Effect\_of\_environmental\_resilience\_on\_processing\_disturbance\_events = GRAPH (Environment\_stabilization\_ratio) Points: (0.000, 0.000), (0.100, 0.360), (0.200, 0.660), (0.330, 1.000), (0.400, 1.100), (0.500, 1.200), (0.600, 1.300), (0.700, 1.400), (0.800, 1.600), (0.900, 1.800), (1.000, 2.000) UNITS: dmnl Effect\_of\_individual\_resilience\_on\_processing\_disturbances = GRAPH(Psychological \_Resilience) Points: (0.000, 0.0133857018486), (0.100, 0.0359724199242), (0.200, 0.0948517463551), (0.300, 0.238405844044), (0.400, 0.53788284274), (0.423, 1.000), (0.600, 1.46211715726), (0.700, 1.76159415596), (0.800, 1.90514825364), (0.900, 1.96402758008), (1.000, 1.98661429815) UNITS: dmnl Effect\_of\_temporary\_disturbance\_on\_increasing\_disturbance\_events = GRAPH(Disturbance\_to\_temporary\_adjustments\_ratio) Points: (0.000, 0.0133857018486), (0.400, 0.0359724199242), (0.800, 0.0948517463551), (1.200, 0.238405844044), (1.600, 0.53788284274), (2.000, 1.000), (2.400, 1.46211715726), (2.800, 1.76159415596), (3.200, 1.90514825364), (3.600, 1.96402758008), (4.000, 1.98661429815) UNITS: dmnl Environment\_stabilization\_ratio = (Environment\_Stabilized\_Adjustments)//(Environment\_Stabilized\_Adjustments+Environment\_Temporary\_Adjustments) UNITS: dmnl Increasing\_disturbance\_events = Effect\_of\_temporary\_disturbance\_on\_increasing\_ disturbance\_events\*Number\_of\_disturbance\_events\_per\_month UNITS: event/month

Initial\_disturbance\_events = Normal\_disturbance\_events\_per\_month\*Time\_needed \_to\_resolve\_events\*2

```
UNITS: event
    Initial_Disturbance_ratio_input = INT(Environment_Temporary_Adjustments//Max
_adjustment_depending_on_resiliency)
        UNITS: dmnl
    Max_adjustment_depending_on_resiliency = Disturbance_events_that_can_be_dealt
_with_resiliency*Number_of_adjustments_needed_per_disturbance_events
        UNITS: adjustment
    Multiple_shocks_to_the_system = (1+"MULTIPLE_SHOCKS_Switch_1=_crisis_on"
*SUM(Multiple_Shocks))
        UNITS: dmnl
    Normal_disturbance_events_per_month = 2.5
        UNITS: event/month
    Number_of_adjustments_needed_per_disturbance_events = 3
        UNITS: adjustment/event
    Number_of_disturbance_events_per_month = Normal_disturbance_events_per_month
*Multiple_shocks_to_the_system
        UNITS: event/month
    Time_needed_to_resolve_events = 2
        UNITS: month
    Total_adjustments_needed = (Disturbance_events*Number_of_adjustments_needed
_per_disturbance_events)
        UNITS: adjustment
    **********
    Policy_switches:
    **********
    CAPACITY_RESPOND_STWITCH = 1
        UNITS: dmnl
    Intensity_of_P1_on_organizational_response_time = 3
        UNITS: dmnl
    Intensity_of_P1_on_revisiting_adjustments = 2
        UNITS: dmnl
    Intensity_of_P1_on_temporary_adjustment_time = 2
        UNITS: dmnl
    Intensity_of_P2_on_fraction_of_dropping_out = 2
        UNITS: dmnl
    Intensity_of_P2_on_fraction_of_using_services = 1.6
        UNITS: dmnl
    Intensity_of_P2_on_time_to_hire = 4
        UNITS: dmnl
    Intensity_of_P2_on_waiting_time = 2
        UNITS: dmnl
    Intensity_of_P3_on_distress_tolerance = 4
        UNITS: dmnl
    Intensity_of_P3_on_stabilised_adjustment = 4
        UNITS: dmnl
    P1:_Environment_Fast_Adaptation = 0
        UNITS: dmnl
    P2_Mental_Wellbeing_Service_Fast_Response = 0
        UNITS: dmnl
    P3_Collective_growth = 1
        UNITS: dmnl
    **********
    Psychological_resilience_structure:
```
\*\*\*\*\*\*\*\*\*\*

Accumulative\_people\_stressed\_up(t) = Accumulative\_people\_stressed\_up(t - dt) + (stressing\_up\_flow) ∗ dt

INIT Accumulative\_people\_stressed\_up = 0

UNITS: people

Disturbance\_to\_resilience\_ratio = Disturbance\_events//Disturbance\_events\_that\_can \_be\_dealt\_with\_resiliency

UNITS: dmnl

Dropping\_out\_services = IF P2\_Mental\_Wellbeing\_Service\_Fast\_Response=1 THEN "Highly\_stressed\_people\_with\_HW\_services"\*Fraction\_of\_dropping\_out\_per\_month/ Intensity\_of\_P2\_on\_fraction\_of\_dropping\_out ELSE Fraction\_of\_dropping\_out\_per\_month \*"Highly\_stressed\_people\_with\_HW\_services"

UNITS: people/month

Fraction\_of\_dropping\_out\_per\_month = 0

UNITS: 1/month

Fraction\_of\_self\_recover = 0.2

UNITS: dmnl

Fraction\_of\_service\_using\_among\_highly\_stressed\_people = IF P2\_Mental\_Wellbeing \_Service\_Fast\_Response=1 THEN Normal\_fraction\_of\_service\_using\_among\_highly \_stressed\_people\*Intensity\_of\_P2\_on\_fraction\_of\_using\_services ELSE Normal\_fraction \_of\_service\_using\_among\_highly\_stressed\_people

UNITS: dmnl

Fractional\_effect\_of\_tolerating\_emotional\_distress\_through\_recovering\_process\_per \_month = ("Self-\_recovering"+Recovering\_by\_using\_services)\*Fractional\_emotional\_ tolerance\_acquired\_from\_recovering//Sum\_of\_population\_stocks

UNITS: 1/month

Frequency\_of\_sessions\_attended\_per\_month\_per\_people = 4

UNITS: session/people/month

Highly\_stressed\_people\_reaching\_out\_services = "Highly\_Stressed\_people\_without \_HW\_services" \*Fraction\_of\_service\_using\_among\_highly\_stressed\_people

UNITS: people

Highly\_stressed\_people\_that\_can\_be\_scheduled = Remaining\_monthly\_sessions\_ capacity//Frequency\_of\_sessions\_attended\_per\_month\_per\_people

UNITS: people

"Highly\_Stressed\_people\_with\_HW\_services"(t) = "Highly\_Stressed\_people\_with \_HW\_services"(t - dt) + (Using\_services - Recovering\_by\_using\_services - Dropping\_out \_services) ∗ dt

INIT "Highly\_Stressed\_people\_with\_HW\_services" = INIT\_Highly\_stressed\_ using\_HW\_services

UNITS: people

"Highly\_Stressed\_people\_without\_HW\_services"(t) = "Highly\_Stressed\_people\_ without\_HW\_services"(t - dt) + (Stressing\_up + Dropping\_out\_services - "Self-\_recovering" - Using\_services) ∗ dt

INIT "Highly\_Stressed\_people\_without\_HW\_services" = INIT\_highly\_stressed \_not\_using\_HW\_services

UNITS: people

"Low/mild\_stress\_people"(t) = "Low/mild\_stress\_people"(t - dt) + ("Self-\_recovering" + Recovering\_by\_using\_services - Stressing\_up) ∗ dt

INIT "Low/mild\_stress\_people" = INIT\_low\_stress\_population

UNITS: people

Normal\_fraction\_of\_highly\_stressed\_symptoms = 0.2

UNITS: dmnl

Normal\_fraction\_of\_service\_using\_among\_highly\_stressed\_people = 0.5 UNITS: dmnl

Recovering\_by\_using\_services = "Highly\_Stressed\_people\_with\_HW\_services"// Time\_to\_recover UNITS: people/month "Self-\_recovering" = "Highly\_Stressed\_people\_without\_HW\_services"\*Fraction\_of \_self\_recover//"Time\_to\_self-recover" UNITS: people/month Stressing\_up = "Low/mild\_stress\_people"\*Disturbance\_to\_resilience\_ratio\*Normal \_fraction\_of\_highly\_stressed\_symptoms//Time\_to\_change\_stress\_level UNITS: people/month Stressing\_up\_flow = Stressing\_up UNITS: people/month Time\_to\_recover = Total\_number\_of\_service\_sessions\_needed\_per\_person//Frequency\_of\_sessions\_attended\_per\_month\_per\_people UNITS: month "Time\_to\_self-recover" = 3 UNITS: month Time\_to\_change\_stress\_level = 6 UNITS: month Using\_services = MIN(Highly\_stressed\_people\_that\_can\_be\_scheduled, Highly\_ stressed\_people\_reaching\_out\_services)/Target\_Waiting\_Time

### UNITS: people/month

#### **References**

## *Article* **Addressing Parameter Uncertainty in a Health Policy Simulation Model Using Monte Carlo Sensitivity Methods**

**Wayne Wakeland 1,\* and Jack Homer <sup>2</sup>**

**Abstract:** We present a practical guide and step-by-step flowchart for establishing uncertainty intervals for key model outcomes in a simulation model in the face of uncertain parameters. The process starts with Powell optimization to find a set of uncertain parameters (the optimum parameter set or OPS) that minimizes the model fitness error relative to historical data. Optimization also helps in refinement of parameter uncertainty ranges. Next, traditional Monte Carlo (TMC) randomization or Markov Chain Monte Carlo (MCMC) is used to create a sample of parameter sets that fit the reference behavior data nearly as well as the OPS. Under the TMC method, the entire parameter space is explored broadly with a large number of runs, and the results are sorted for selection of qualifying parameter sets (QPS) to ensure good fit and parameter distributions that are centrally located within the uncertainty ranges. In addition, the QPS outputs are graphed as sensitivity graphs or box-and-whisker plots for comparison with the historical data. Finally, alternative policies and scenarios are run against the OPS and all QPS, and uncertainty intervals are found for projected model outcomes. We illustrate the full parameter uncertainty approach with a (previously published) system dynamics model of the U.S. opioid epidemic, and demonstrate how it can enrich policy modeling results.

**Keywords:** simulation model; uncertainty analysis; optimization; sensitivity testing; Monte Carlo randomization; opioid epidemic

### **1. Introduction**

#### *1.1. Background and Approach*

System dynamics (SD) models frequently employ parameters for which solid empirical data are not available. Modelers often use expert judgment to provide estimates of parameter values for which empirical data are not available. When several experts are available, and a formal process (such as Delphi) produces convergent estimates, modelers may have reasonable confidence in the parameter values despite the lack of data. Nevertheless, these parameter values remain uncertain to a degree—as is indeed true even for measured parameters, due to issues including small sample sizes and definitional variation.

To address parameter uncertainty, SD modelers test alternative parameter values in order to understand their degree of influence in the model. Formal sensitivity analyses can be run using features provided in popular SD software packages, and results can be displayed in a table or portrayed graphically as a "Tornado diagram" [1].

Once the modeler has identified influential parameters, additional effort is applied, as time and budget allow, to increase confidence that the values of these influential parameters are well supported. When reliable empirical data are available, ideally from multiple sources, the parameter value may be fixed and used with confidence. Usually, however, some, perhaps even many, parameter values typically remain uncertain.

This paper describes how the model analysis features available in the VensimTM software package (as well as other popular packages including Stella®), can be used to

**Citation:** Wakeland, W.; Homer, J. Addressing Parameter Uncertainty in a Health Policy Simulation Model Using Monte Carlo Sensitivity Methods. *Systems* **2022**, *10*, 225. https://doi.org/10.3390/ systems10060225

Academic Editor: Peter Jones

Received: 17 September 2022 Accepted: 12 November 2022 Published: 18 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

represent and incorporate parameter uncertainty into the analysis of SD model results, including providing outcome uncertainty intervals.

In this paper, we present a flowchart describing a step-by-step process for incorporating parameter uncertainty in SD models. We then illustrate use of the method using a previously published model of the U.S. opioid epidemic [2]. We close with a discussion of how model results under uncertainty can be described for interested parties such as policy makers or other end users of the analyses.

### *1.2. Prior Work in This Area*

Ford [3] and Sterman [4] discussed uncertainty analysis vis-à-vis dynamic models. Helton et al. [5] provided a general overview of sampling-based methods for sensitivity analysis, including traditional Monte Carlo (TMC). Dogan [6,7] discussed confidence interval estimation in SD models using bootstrapping and the likelihood ratio method. Cheung and Beck [8] explained Bayesian updating in Monte Carlo processes, and background material on the mathematics of Markov Chain Monte Carlo (MCMC) is readily available [9,10].

Background on sensitivity analysis methods for dynamic modeling, and SD in particular, including the TMC and MCMC methods being employed in this paper, may be found in Fiddaman and Yeager [11], Osgood [12], Osgood and Liu [13], and Andrade and Duggan [14]. Using TMC methods to search a parameter space is sometimes considered to be a brute force or totally random search, whereas the MCMC method uses a Bayesian update process to guide the search of parameter space. In at least some contexts, it has been shown that MCMC optimally creates statistically valid samples [9].

Although publications describing the application of these methods are plentiful in some scientific and engineering disciplines, publications featuring their use with SD models are scant. A search for "system dynamics" AND ("MCMC" or "monte carlo") returned only a handful of publications. Five relevant examples are Jeon and Chin [15], who described their use of TMC with an SD model of renewable energy; Sterman et al. [16], who applied MCMC in a model of bioenergy; Garfazadegan and Rahmandad [17] who used MCMC to estimate parameter values for a COVID-19 model, Lim et al. [18] who applied MCMC in a model of the U.S. opioid crisis; and Rahmandad et al. [19] who applied MCMC in another model of COVID-19.

#### **2. Materials and Methods**

#### *2.1. The Process: Initial Steps*

Figure 1 presents the initial steps of a process for incorporating uncertainty analysis into SD models, with further steps shown in Figures 2–4. The complete process is shown in a single flowchart in Supplement Part S1.

The process starts at *Create model* **&** *modify as needed*. This would be a model that employs uncertain parameters and includes dynamic outcome variables that strive to match the dynamics seen in real world reference behavior data. To use the methods described in Figure 1, one needs to *Define error metric variables* and *Add error metrics to model*. A useful example could be to compute the mean absolute error (MAE) between the model calculated time series and the reference behavior time series for each outcome variable. Care must be taken to consider how to compute MAE when reference data are incomplete so as not to distort results. When different outcomes have very different scales, it is useful to use MAEM, which stands for MAE over M, where M is the mean value of the metric. In addition, composite error statistics are added to the model, such as the average of the MAEMs over all the outcomes, and the maximum value of the individual MEAMs. These are used later in the process for identifying well-fitting parameter sets. There are statistical macros available for Vensim to help with this (Supplement Part S2).

**Figure 1.** Process for addressing simulation model uncertainty, initial steps.

Another initial step is to *Estimate parameter uncertainty ranges*, a lower and upper bound for each uncertain parameter, and to *Specify weights for the outcome variables*. These are needed for the algorithm used in a key part of the process: *Estimate uncertain parameter values*. This process employs a Powell optimization process that uses an objective function consisting of model outcomes vs. reference data. For more information, see Menzies et al. [20] for a tutorial on using Bayesian methods to calibrate health policy models. Each term of the objective function represents one of the outcome time series of interest. The algorithm strives to minimize the differences between model and reference data. Weights are needed because outcomes may have very different scales and variances. This tends to be an iterative process, so Figure 1 contains a feedback loop, and some of the connections are bidirectional. The end result of this step is a set of optimized parameter values for the uncertain parameters (OPS). Its average MAEM might be 0.1 and its maximum MAEM (worst MAEM for any one of the outcomes) might be 0.2.

#### *2.2. The Process: Intermediate Steps for the Traditional Monte Carlo (TMC) Approach*

At this point, the user may elect to use traditional Monte Carlo (TMC) or Markov Chain Monte Carlo (MCMC). TMC is discussed first; see Figure 2.

*Make very large traditional Monte Carlo (TMC) run* employs Vensim's sensitivity feature to perform a very large number of model runs, millions if there are many uncertain parameters. Using this feature requires the user to specify how many runs to perform, the seed to start with, the type of sampling (e.g., multivariate, Latin hypercube, etc.), and which parameters to vary and how. We used multi-variate, which is a totally random search process. Latin hypercube strives to cover a large parameter space more efficiently. Both may be valid choices. We experimented with Uniform but settled on Triangular with the mode specified as the value from the OPS. Since we were focused on overall error, we changed the output save period (SAVEPER) to be the length of the run. This kept the output file size manageable. We also specified some additional variables to be saved (all the parameters being varied are automatically saved). These additions were the average MAEM and maximum MAEM for the run.

**Figure 2.** Intermediate steps for the traditional Monte Carlo (TMC) approach.

Note that when more than one million runs are needed, we find it practical to perform one million at a time and to change the seed for each run. At the end of a sensitivity run, the data output file (.vdf in Vensim) is exported to a tab-delimited file for further analysis.

Next, the user needs to *Determine cutoffs for a qualified parameter set (QPS)* which involves determining how well-fitting a candidate parameter set must be to warrant its inclusion in the QPS. The vast majority of the TMC runs will not be well-fitting in the manner of the OPS. One rational approach would be to accept parameter sets that performed nearly as well as the OPS. For example, for a model in which the OPS has an average MAEM = 0.1 and the maximum MAEM = 0.2, the cutoffs for the inclusion of a candidate parameter set could be average MAEM < 0.12 and maximum MAEM < 0.25.

*Sort TMC runs* **&** *create N-QPS* employs Excel's data/import external data from file to read in the TMC results. The result will be M columns by N rows, where M is the number parameters being varied + K, and N is the number of runs. K is the number of saved variables times the number of saved times per variable, which could be one, two (start time and end time), or more. The results are then sorted by average MAEM and all rows > average MAEM cutoff are discarded. The remainder is sorted by maximum MAEM and all rows > maximum MAEM are discarded. This will likely leave a very tiny fraction of runs, perhaps a few hundred out of a million. These N runs are the qualified parameter sets (N-QPS).

#### *2.3. The Process: Intermediate Steps for Markov Chain Monte Carlo (MCMC) Approach*

A very different approach for creating a set of runs to be used to evaluate the impacts of parameter uncertainty is to use the MCMC approach; see Figure 3.

This begins with *Estimate Optimal Outcome Weights*. First, one adds to the model a weight variable for each outcome, instead of specifying the weights numerically as in the MC process. One also needs to specify the search range for outcome weight constants. The search ranges for the uncertain parameters can be the same as for the MC approach. Or, these ranges could be broader than those used with the prior approach. No sampling distribution such as Uniform or Triangular is specified because the search of parameter space is guided by a heuristic, not by random, sampling. Guidance regarding how broad to set the ranges for MCMC varies by expert. We have heard from some, the comment that very broad ranges will give the algorithm "room to work". Others have suggested that ranges should not include implausible values. Both comments are sensible, suggesting that a careful study of this question is needed. To proceed, first a Powell search is performed

that optimizes the uncertain parameters and the weights simultaneously. In our experience, allowing larger ranges may allow the algorithm to find an optimum that achieves slightly lower average MAEM, but higher maximum MAEM. The reason is that some of the individual outcomes are sacrificed to achieve lower overall error. However, these results may be less realistic.

**Figure 3.** Intermediate steps for the Markov Chain Monte Carlo (MCMC) approach.

Next is to *Conduct MCMC search* **&** *generate SVS*. MCMC, in the manner of MC, varies the uncertain parameters over the specified ranges, but uses in the objective function the outcome weights that were found in previous step. The Markov Chain-driven search algorithm is designed to create a statistically valid sample (abbreviated here as SVS). It can be quite large and contain many duplicates (which apparently helps to assure a sample with the correct properties). Model fitness may be further improved, but be less plausible, as mentioned earlier. Results are then exported as a tab-delimited file.

One then must *Import results into Excel*, and, if the population of runs produced by the MCMC search is excessively large, one can select a random sample of desired size M from the SVS. This subsample of size M, which is still a statistically valid sample (M-SVS), is comparable to the N-QPS set of runs created by the conventional MC process described earlier. One can proceed to computing statistics by parameter using the M-SVS and/or use it to run file-driven sensitivity runs, analyze alternative scenarios, etc.

#### *2.4. The Process: Final Steps for Both TMC and MCMC*

The final steps of the process are shown in Figure 4.

*Save sensitivity (TMC or MCMC) parameters* as a tab-delimited txt file creates the file needed for file-driven sensitivity runs.

Within Excel, the user can proceed to *Compute/graph statistics by parameter* to see what the values in N-QPS (or M-SVS) file are for each parameter, for example, to see if the entire range of possible values for a given parameter is represented in N-QPS, if the mean of this sample is near the value of the parameter in the OPS. Additionally, it could be used to determine the confidence interval of the estimate for the mean of the parameter based on this sample and examine the distribution (shape) of the sample. Does this information raise any red flags with respect to the OPS or N-QPS? Or, does this information indicate that the N-QPS may be representative of the entire parameter space?

**Figure 4.** Final steps of the process.

Next, the N-QPS (or M-SVS) file can be used in two primary ways. One is the *Run file-driven Trajectory Sensitivity Run for Baseline*, for which the user can change the SAVEPER to a useful time-period, such as by year, in order to create a series of N (or M) trajectories for each outcome. The result is exported to a tab-delimited file that can be imported into Excel to *Compute/graph baseline statistics by primary outcome*, including uncertainty intervals for the model calculated outcome at each time period, since there is now a sample of N model-calculated values at each time point. Excel can also be used to create box and whisker plots at each time point, with or without outliers. However, we found it more expedient to create these plots using Python.

For the other primary use of the N-QPS (or M-SVS) file, explained below, the user first needs to *Specify primary high-level outcomes of interest* such as total cost and net performance. The previous sensitivity run was focused on behavior over time for all outcome time series, but to compare alternative scenarios in an overall sense, a few key end-of-run metrics are needed. In addition, one also needs to *Specify alternative scenarios of interest*; these might be different configurations or policy changes that could be implemented in the model as switches that are used in conjunction with magnitude-of-impact parameters and timing parameters linked to specific model constants.

Using the N-QPS or M-SVS, the user will change the SAVPER back to the End of Run and then *Make file driven sensitivity runs for baseline and each alternative scenario*. The results, for the baseline and for each scenario, are matched pairs of data points (where all of the uncertain parameter values are the same for baseline and the alternative). This means that the distributions of the differences can be used to determine, in a statistically valid fashion, the credible interval estimates (a term from Bayesian statistics), for the mean of the difference by outcome between baseline and the alternative, that can be attributed to parameter uncertainty.

We suggest that one *Assemble a consolidated spreadsheet* to perform these analyses. In the upper left corner, is O × N, where O is the number of overall outcomes and N is the number of runs. Note that the raw sensitivity results file has columns for all the parameters as well as columns for each outcome. The user selects and copies only the end columns for each outcome into the consolidated sheet. The data for Alternative 1 is placed below the Baseline results; Alternative 2 below Alternative 1, etc. Then, starting at the first row of Alternative 1, columns will be added to compute the differences between the Alternative 1 numbers and baseline number, one column for each outcome. Similarly for Alternative 2 vs. baseline, and so forth. Columns for percentage differences can also be added.

Finally, *Create summary table(s) or tab(s)* is used to provide the results of relevant calculations, such as means and their credible intervals for each outcome at baseline and for each alternative. In addition, the results for the differences by outcome between baseline and Alternative 1, baseline and Alternative 2, etc., may be compiled. Similarly, for the percentage differences.

#### *2.5. About the Opioid Epidemic Model*

To illustrate the process and results, the TMC method was applied to an SD model of significant complexity that explored the opioid epidemic in the United States from 1990 to 2030 including the impacts of alternative policies (hereafter the Opioid Epidemic Model) [2,21]. SD has been frequently applied to drug abuse and other areas of public health and social policy [22–26]. The current model adopts the basic scientific approach and some of the same elements as these forerunner models, but the current model was completely redeveloped to address current needs.

A committee of the National Academies of Sciences/Engineering/Medicine considered the complexities of the opioid epidemic and stated that for informed decision making "a true systems model, not just simple statistics" was needed because "decisions made about complex systems with endogenous feedback can be myopic in the absence of a formal model" [27]. The NAS cited the system dynamics model by Wakeland et al. [24] as an example of a true systems model.

Best practices for model development, testing, and reporting are well documented [4,28–30]. The model building process for the Opioid Epidemic Model involved amassing evidence from many different sources and developing a dynamic structure that reproduces historical trends and is sufficiently parsimonious to allow explanation. Uncertainties are addressed through sensitivity testing and the multivariate TMC testing described previously that allows results to be reported with rigorous credible intervals.

A complete list of the approximately 340 interacting model equations (including about 80 input constants, 17 time series for input and validation, and 240 output variables) is available in the model's reference guide [21]. This guide also provides stock and flow diagrams for each sector of the model and describes how model constants and time series were estimated, including data sources. The model time horizon was 1990–2030. Figure 5 shows the basic model structure.

**Figure 5.** Opioid Epidemic Model diagram showing high level stocks, flows, and primary feedback structures. Red lines are causal links with negative valence, while blue lines have positive valence. Reproduced by permission from Homer and Wakeland 2020 [2].

#### **3. Results**

#### *3.1. Error Metrics*

As the first step in the TMC approach, error metrics were created, in this case, MAEM (mean absolute error over mean) for multiple outcome time series with historical data counterparts. We used the Vensim SSTATS macro provided by Professor John Sterman at MIT; see Supplement Part S2. This code was inserted via text editor (we used Notepad) at the top of the model file starting at the second line. Supplement Part S3 presents the model code for calculating the statistics, which was also inserted into the model file via text editor, typically after the Control block, which is located after the user-defined model variables/constants section. We utilized an Excel workbook that included a worksheet with Time in Row 1 and historical time series data in the rows below. In some cases, there were data only for some of the historical years. The SSTATS macro calculations are designed to handle this correctly. An example Vensim equation for reading one these data time series is GET XLS DATA ('model RBP data.xlsx', 'RBP', '1', 'B2').

Another step to prepare the model for analysis was to add a custom graphs file via the Control Panel and add custom tables that display the variables calculated by the SSTATS macro (R2, MAPE, MAEM, RMSE, Um, Us, Uc, and Count), one table for each Outcome variable. I/O Objects were added to a model view to display the custom tables.

Figure 6 shows a sample of these statistical parameter displays created using the SSTATS macro, with the MAEM statistics on the left and R-squared statistics on the right. These outputs were from the optimized parameter set (OPS) version of the Opioid Epidemic Model; see discussion below. In addition to the weighted average of all the MAEMs of 8.9% (bottom left 2020 column), it is clear that OD death-related trajectories were reproduced more accurately than, for example, heroin-related trajectories. The MAEM for *OD deaths total* was 3.6%, whereas for *H initiates* it was 18%. This large difference was due mostly to noise in the data representing the actual values. OD death data are the actual data collected by the CDC, and the changes from year to year are modest. Data on persons initiating heroin in a given year is based on a small sample of a difficult-to-access population and the resulting estimates by year vary dramatically. In fact, the difference between the data and the corresponding three-period moving average was 12.5%, suggesting that most (70%) of the 18% model-calculated MAEM was due to noise in the data.

**Figure 6.** Selected SSTATS displays (MAEM and R-squared for outcome variables).

#### *3.2. Optimized Parameter Set (OPS)*

Before using optimization to estimate uncertain parameter values, the user must specify an appropriate range for each parameter to be used by the Powell search algorithm. Supplement Part S4 provides an example Vensim optimization control file (.voc) which specified the uncertain parameters to be varied. Table 1 shows the first few rows of the model parameter spreadsheet by parameter, including the minimum and maximum values. Many of the 80 model parameters had empirical support from the literature, which helped to reduce the uncertainty ranges for these parameters.

In addition, weights must be specified for the objective function, informed by the relative magnitudes and variances of the outcome variables. In this example, weights were set so that metrics with solid empirical data were given more weight than those for which empirical data were scant or less reliable.

The optimization runs using the CG (calibration Gaussian) option ran for several minutes on a laptop computer. A blend of modeler judgment and optimization results was used iteratively to select the values for the optimized parameter set. The MAEMs for many of the outcomes were reduced compared with a purely manual calibration process. The weighted (based on the amount of data available for each outcome) average of the MAEMs was 8.9%, with the largest individual MAEM being 17.9%.

**Table 1.** Portion of the opioid model parameter spreadsheet showing parameter name, units, value (central estimate), source (or optimized), and the minimum and maximum values used during optimization runs.

Figure 7 shows the trajectories for two example outcomes calculated by the model using the optimized parameter set compared with the available historical data for these metrics. The MAEMs for the two outcomes in Figure 7 were 9% and 4%. The calculated *Persons with Addiction* trajectory was somewhat biased downward, whereas the *OD Deaths* trajectory matched very well. One possible reason could be the amount of data available: the standardized errors for the first variable may carry somewhat less weight than the second one. Additionally, note that measurement uncertainty was greater for the first variable than the second one.

**Figure 7.** Model fit to two key outcomes using the optimized parameter set.

### *3.3. Qualified Parameter Sets (N-QPS)*

Creating a large Monte Carlo run is easy using the sensitivity feature in Vensim, which uses or creates a .vsc file. One option for exploring parameter space could be to use Uniform distributions between thoughtfully chosen minima and maxima. These values could be the same as those used during calibration, or perhaps narrowed somewhat. Another option could be to use a Triangular distribution with the optimum value as the mode, which would increase samples near the optima. For this illustration, Triangular distributions were used. The user also needs to specify a .lst file that lists which variables are to be saved for each run. Sensitivity runs automatically save the varied parameters used for each run, and in addition, it is useful to know the maximum MAEM and average MAEM (to be used to select which runs are qualified). Since we do not need to know the time trajectories, we set the model SAVEPER to 40. This kept the output file modest in size. Supplement Part S5 shows the .vsc and .lst files.

The simulation was set to create one million runs, which took about seven hours on a laptop. Ten such runs were made, changing the random seed for each run, creating a total sample of ten million runs. Each resulting data output file (.vdf) was exported (via the model menu) to a tab-delimited file and read into Excel using data/get external data/from text (switch to All Files to see .tab files). This was sorted by weighted average of MAEMs. The weighted average MAEM was below 0.11 for 300–600 runs in each of the one million runs, and these rows were kept. The other rows were deleted.

The file was then sorted by maximum MAEM, and 100–130 runs were less than 0.20 and kept. The rest were deleted. These ten spreadsheets were combined, yielding a sample of 1119 runs that accomplished <11% weighted average error and <20% maximum error. Rows were added below the sample to calculate useful statistics about the sample, regarding both the distributions of the input parameters and the how well each run performed. Table 2 shows the first and last few rows of the spreadsheet.

**Table 2.** First and last few rows and columns of the file used to create the N-QPS.

We next examined the properties of the parameter samples contained in the N-QPS. The rows at the bottom of Table 2 show the minimum, maximum, mean, and standard deviation of each parameter. For more detail, histogram plots were created, with examples shown in Figure 8. The parameter shown on the left, *Addicted PONHA move to heroin rate initial*, was limited to be in the 0.01 to 0.03 range. As the histogram shows, nearly all of

this allowed range was included in the N-QPS. The mode was about 0.0235, which was slightly higher than its value of 0.021 in the optimized parameter set. The parameter on the right, *Consumption mgs ME per addicted PO abuser per month initial*, was limited to be in the range of 2000 to 5000; the values in the N-QPS were well inside these limits (2092 to 4566). The mode was 3330, very slightly higher than its value of 3200 in the optimized parameter set. One can examine all of the parameter histograms in this fashion to build confidence in the sample.

**Figure 8.** Histograms for two uncertain parameters resulting from the TMC approach.

For comparison, Figure 9 shows a histogram that results from using the MCMC algorithm rather than TMC, again considering the parameter that was on the left side of Figure 8. For MCMC, nearly all of the samples fell in the range from 0.02 to 0.0216 with a width of 0.0016, whereas the range for same parameter in Figure 8 is from 0.014 to 0.03 with a width of 0.026 (greater by 16x), and the shape has a sharp rather than rounded peak. We found that nearly all of the parameter distributions from testing MCMC on the Opioid Epidemic Model had this narrow and pointed characteristic. The problem was not skewness, but the implausibly high degree of precision in the estimated parameter value, within just a few percentage points.

Using different settings for the MCMC algorithm yielded similar very narrow and pointed samples, prompting the decision to rely on the TMC-based method that tests a much broader range of potential values for these parameters. However, other researchers have reported successful use of the MCMC method [15–19], indicating that the MCMC approach may be preferred in many cases.

#### *3.4. Sensitivity Runs for the Baseline Scenario*

To put the N-QPS sample from the MC-based method to use, file-driven sensitivity runs were performed. The SAVEPER was changed to 40, and the Sensitivity tool was selected. Previous settings were cleared, and the Type of sensitivity run was changed to File. The N-QPS text file was selected as the file. Next a save list was created. The previous save list (with the statistical variables) was cleared, and three primary outcome variables were entered: *Total opioid addicts, Opioid overdoses seen at ED, and Total opioid OD deaths.* This run took just a minute to complete, and the data output file (.vdf) it produced was exported via Vensim and imported into Excel.

**Figure 9.** Histogram for an uncertain parameter resulting from the MCMC approach. Note that all of the bars in this diagram are contained within the fourth bar on the left side of Figure 8.

The first 50 columns of this trajectory spreadsheet provided the run number (1–1119) and the values of the 49 uncertain parameter values used for each specific run. The remaining 41 columns were values from time 1990 to 2030 for the three outcome variables. Figure 10 shows the trajectory under uncertainty for *Overdoses Seen at ED* using a boxand-whiskers plot. One can also use the data in the spreadsheet to report useful sample statistics, such as mean and standard deviation.

**Figure 10.** Excel box-and-whisker time series plot for *Overdoses Seen at ED* under the baseline scenario.

Unfortunately, there was no easy way using Excel to superimpose the historical data on this plot. Python software could be used to create such a plot as needed. Figure 11 shows the plot produced using Python for *Total opioid OD deaths*, showing the credible interval trajectory overlaid with green dots for the historical data. The credible interval was similar in spirit to the confidence intervals used to characterize statistical samples from actual populations.

**Figure 11.** Python box-and-whisker time series plot for *Total Opioid OD Deaths* under the baseline scenario, overlaid with historical data (green dots).

To prepare the file needed by Python, a version of the trajectory spreadsheet was created, named Outcome Trajectory Data for Python Plotting (OTDPP), in which the columns for run number and the parameter values were deleted. Additionally, two rows were added below the first row. Row 2 should contain the actual data to be overlaid. For each outcome variable, the appropriate row of the RBP spreadsheet was copied above its corresponding section in OTDPP. Since the RBP data spreadsheet contained up to 30 values of historical data for each outcome trajectory, and since the data sections for each outcome in OTDPP were oriented horizontally, it was a simple copying process, even for several outcomes. Row 3 in OTDPP should provide the time values, which was easily created in Excel or could be copied from the historical data spreadsheet (41 cells copied once and pasted N times in the OTDPP, where N is the number of outcome variable, 3 for this example).

#### *3.5. Sensitivity Runs for Alternative Scenarios (Policy Testing under Uncertainty)*

Another use of the N-QPS is to compare a baseline case with policy runs, to determine how much impact parameter uncertainty may have with respect to the most important outcome indicators. In the case of the Opioid Epidemic Model, three key indicators were the number of people with opioid addiction (also known as opioid use disorder or OUD), opioid overdose events, and opioid overdose deaths. These metrics are saved for each of the parameter sets in the N-QPS, for the baseline condition and for each alternative. The net change in outcome, and policy vs. baseline, were calculated run by run, in both absolute and percentage terms.

The Opioid Epidemic Model was run using the file-driven sensitivity tool reading in the 1119-QPS parameter sets, and for five potential policies aside from the baseline. The final values for the three key metrics were saved for each of the 1119 runs, and a

spreadsheet was created to compute the differences in outcome metrics, run by run. This spreadsheet had 1119 rows and 48 columns. The first three columns were the final baseline values for each of the three metrics. Then, nine columns were created for each of the five policies. These were the three columns of raw outcomes, plus three columns that calculated the change, and three more that calculated the percentage change. A row was added at the bottom to compute the mean percentage change in each outcome due to each policy. Two more rows computed minimum change in each metric and the maximum change in each metric.

Table 3 presents a summary of the parameter uncertainty analysis of the impacts of the various policies. The *% change vs. Baseline* under OPS may be compared with the *Mean* %∆ and the range (min, max) of percentage changes under the QPS 1119. Even though the mean of the QPS-driven change was generally close to the results using the OPS, seeing the uncertainty interval provides more information to inform policy decisions. Three examples are highlighted in yellow, as follows:

**Table 3.** Impact of parameter uncertainty on the simulated results of policy changes.

#### **4. Discussion**

Here we have demonstrated a practical method for analyzing the effects of parameter uncertainty on simulated model projections, using a health policy example. The process summarized in Figures 1–4 (and Supplement Part S1) provides two alternative approaches: traditional Monte Carlo (TMC) and Markov Chain Monte Carlo (MCMC). We featured the former in this paper, because of its greater familiarity and ease of understanding.

However, it should be noted that, despite our difficulties with MCMC (specifically, the overly narrow parameter distributions we observed), it is the theoretically superior and less time-consuming of the two approaches. Tom Fiddaman at Ventana Systems (the makers of Vensim) has suggested that a better choice of likelihood function for MCMC might have produced better dispersed parameter distributions.

Indeed, other research teams have had success with MCMC [15–19]. Although MCMC yields theoretically superior results, further examination and characterization of the parameter distributions resulting from MCMC seems to be in order. Garfazadegan and Rahmandad [17] note in their Appendix A "rather tight confidence intervals coming from MCMC methods directly applied to large nonlinear models," which they address via heuristics that scale the likelihood function.

A possible limitation of both TMC and MCMC is their reliance on the goodness of fit to historical data. Forrester [32] warned that "the particular curves of past history are only a special case." The implication is that this method should not be used to make claims of precision in predicted outcomes, but rather to better appreciate the range of possible outcomes and, especially, the uncertainty in projected policy impacts.

#### **5. Conclusions**

We have presented a step-by-step approach to assessing the degree of uncertainty in simulation model outcomes related to uncertainty in the model's input parameters. Our flowchart summarizes two ways to perform this, and we have demonstrated one of these in detail with a concrete health policy example. Providing uncertainty intervals for the range of possible outcomes from contemplated policy options, as we have described here, could increase the value of simulation models to decision makers.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/systems10060225/s1, Supplement Part S1: Complete process for addressing simulation model uncertainty (uniting Figures S1 to S4) and Supplement Parts S2 to S5.

**Author Contributions:** W.W. and J.H. have met authorship requirements, and contributed to every part of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The original Opioid Epidemic Model was developed under contract in 2019 (see [2]), but the uncertainty analysis described here was performed subsequently and received no funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data presented in this article were extracted from publicly available data sources (see [2]).

**Acknowledgments:** The authors wish to thank Kevin Stoltz who wrote the Python code for graphing uncertainty interval time series, as well as John Sterman and Tom Fiddaman for their advice.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Enabling Mobility: A Simulation Model of the Health Care System for Major Lower-Limb Amputees to Assess the Impact of Digital Prosthetics Services**

**Jefferson K. Rajah <sup>1</sup> , William Chernicoff <sup>2</sup> , Christopher J. Hutchison <sup>3</sup> , Paulo Gonçalves <sup>4</sup> and Birgit Kopainsky 1,\***

**Abstract:** The World Health Organization estimates that 5 to 15% of amputees in any given population have access to a prosthesis. This figure is likely to worsen as the amputee population is expected to double by 2050, straining the limited capacity of prosthetics services. Without proper and timely prosthetic interventions, amputees with major lower-limb loss experience adverse mobility outcomes, including the loss of independence, lowered quality of life, and decreased life expectancy. Presently, the use of digital technology in prosthetics (e.g., 3D imaging, digital processing, and 3D printed sockets) is contended as a viable solution to this problem. This paper uses system dynamics modeling to assess the impact of digital prosthetics service provision. Our simulation model represents the patient-care continuum and digital prosthetics market system, providing a feedback-rich causal theory of how digital prosthetics impacts amputee mobility and the corollary socio-health-economic outcomes over time. With sufficient resources for market formation and capacity expansion for digital prosthetics services, our work suggests an increased proportion of prosthesis usage and improved associated health-economic outcomes. Accordingly, our findings could provide decision support for health policy to better mitigate the accessibility problem and bolster the social impact of prosthesis usage.

**Keywords:** prosthetics; major lower-limb amputations; prosthesis usage; amputee mobility; system dynamics; simulation model; health care system; health policy

### **1. Introduction**

The World Health Organization (WHO) estimates that around 0.5% of any given population require prosthetics and orthotics services [1]. This figure is expected to double by 2050 as a result of ageing populations and rising rates of medical conditions, such as diabetes mellitus, peripheral arterial disease (PAD), and sepsis [1,2]. Particularly for major lower-limb amputations (i.e., above ankle), over 90% of cases in industrialized countries are attributed to PAD (either primary or secondary to diabetes); whereas traumatic injuries make up most cases in developing countries [3–5]. PAD is a progressive vascular disease that commonly causes arterial obstruction in the lower extremities. Known PAD risk factors include cigarette smoking, diabetes mellitus, hypertension, and dyslipidemia, with incidence sharply rising for populations above age 50 [6,7]. PAD progresses to the more severe critical limb ischemia, if not effectively managed at an earlier stage, which could lead to amputation [8].

Major lower-limb amputation, without timely prosthetic intervention, leads to a loss of mobility, which has several ripple effects at both the individual and societal level. It worsens individual health and psychosocial outcomes, including the loss of independence,

**Citation:** Rajah, J.K.; Chernicoff, W.; Hutchison, C.J.; Gonçalves, P.; Kopainsky, B. Enabling Mobility: A Simulation Model of the Health Care System for Major Lower-Limb Amputees to Assess the Impact of Digital Prosthetics Services. *Systems* **2023**, *11*, 22. https://doi.org/ 10.3390/systems11010022

Academic Editors: Jack Homer and Gary B. Hirsch

Received: 13 November 2022 Revised: 27 December 2022 Accepted: 29 December 2022 Published: 3 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

increased depression and self-esteem issues, lowered quality of life, and increased risk of comorbidities and mortality [9–11]. There is also high economic burden on patients, families (increased caregiving), health and welfare systems, as well as the workforce (lower rates of return to work) [1,12,13]. Such negative externalities can be alleviated with the use of prostheses to regain mobility and functional independence [14]. However, WHO estimates that only 5 to 15% of amputees have access to prosthetics services [1]. Even then, approximately 50% of amputees in prosthetics care abandon the process or their prosthetic devices [15,16]. Barriers include high financial costs for treatment, poor health care coverage, prosthetics service capacity constraints, limits in prosthetics technology and fitting, lack of proximity to services, and inadequate continuity of care [14–17].

The recent introduction of digital technology in prosthetics is seen as a viable solution to the accessibility problem [16,18,19]. Digital solutions to prosthesis fitting (henceforth, digital prosthetics) involve a streamlined process of scanning the limb and using a digital software to create a model of the socket for three-dimensional printing. Digital technology reduces the delays, patient time and travel burden, and labor involved in traditional prosthetics. Using traditional methods, the prosthetist must handcraft the socket using plaster casts and test the fittings several times before a definitive socket is manufactured and assembled [19]. With digital prosthetics, this manufacturing delay can be more than halved. In turn, this reduces the chances of the patient's limb and/or weight having changed before receiving the prosthesis device—the main cause of discomfort and pain [16]. The fit challenges are a primary cause for the 50% abandonment [15,16]. Hence, digital prosthetics could lead to higher success rates since the digital design is more accurate, precise, and enables direct translation of a prosthetist's skill-level over a minimum baseline in place; has a much shorter timeframe such that there is little time for limb changes; and results in a more comfortable fit for patients [16].

Moreover, digital prosthetics could improve accessibility by expanding service capacity. With a more streamlined and effective fitting process, each prosthetist can fit more patients in their schedule than it otherwise would have been possible with conventional techniques and processes. Digital technology also frees the prosthetist from their clinic and gives them the flexibility to bring the service to patients through distributed care networks [20]. Accordingly, proponents of digital prosthetics anticipate several positive externalities for amputees, their families, and the economy more broadly. This paper seeks to assess this impact of digital prosthetics service provision on total amputee mobility. Mobility, here, is measured by the proportion of medically eligible amputees who are fitted with a prosthesis and have regained functional mobility. The benefit of digital prosthetics can be further measured by the health-related socio-economic consequences of such mobility; namely, the surplus economic productivity from returning to work and the net economic costs incurred or avoided (health care, family opportunity cost, social and welfare payments).

The purpose of this study is to explore how the adoption of digital prosthetics impact amputee mobility and associated outcomes over time. To assess such changes, we model the key causal mechanisms found in the health care system, including the patient-care continuum and prosthetic service provision. For this purpose, we build and analyze a dynamic simulation model to identify high-leverage points that can enhance the effects of digital prosthetics service provision on mobility outcomes. This paper describes the structure, empirical foundation, illustrative results, and strategic insights from the prototype prosthetic service provision model, as well as how it might be further refined. The results reported in this paper result from two activities: First, we use an in-depth review of the existing knowledge (from literature and expert interviews) coupled with causal loop diagramming [21] to identify core feedback mechanisms driving prosthetic service provision. Second, we developed a formal system dynamics simulation model to characterize the range of outcomes that these processes generate, even in a data-poor context. The end result is an internally consistent theory that provides insights into the determinants of success and failure of digital prosthetics service provision.

#### **2. Materials and Methods**

This study employs the system dynamics (SD) method, using compartment models, for conducting model-based hypothesis testing. SD models seek to simulate and explain problem behaviors by modeling the underlying system structure [21]. Importantly, they offer an "endogenous or feedback perspective" to structural problems [22] (p. 1) that can aid "theory building, policy analysis, and strategic decision support" [23] (p. 11). This endogenous perspective relates to two fundamental methodological tenets: (1) problem behaviors arise from the complex interaction of interrelated components within a closed boundary of a system, and (2) the system components are connected in feedback loops (circular chains of causal relationships), which endogenously generate the observed system behavior [24,25]. In this sense, the SD method "helps construct a causal-loop theory of system behavior in terms of feedback linkages" [26] (p. 400).

SD modeling is well-suited to domains in public health and medicine, with over 300 applications to date – for a review, see [27,28]. The "dynamic complexity in public health" (particularly due to nonlinear effects of multiple interacting variables within the system that affect health outcomes) makes it "difficult to know how, where, and when to intervene" [29] (p. 452). SD simulation modelling can effectively address this challenge and "elucidate the counterintuitive behavior of complex healthcare problems" [28] (p. 1). Particularly for prosthetics provision and related health care policy, SD modelling can support decision-making under uncertainty. The domain of prosthetics services is mired by the lack of robust data collection, contributing to a high level of uncertainty surrounding policy planning [1,16]. SD models, however, can "admit more variables on the basis of logic or expert opinion and for which solid statistical estimates may not be available" [29] (p. 453) and still generate useful insights under such uncertainty.

#### *2.1. Literature Review*

To our knowledge, apart from the two preliminary versions of this work [30,31], there has been no other application of SD to prosthetic service provision or major lower limb amputations in the academic literature. However, this work builds on existing SD literature on emerging medical technologies and innovation diffusion more generally. Paich, Peck and Valant present a model on pharmaceutical product strategy that integrates patient flows, product diffusion and adoption by physicians, and treatment attractiveness to patients [32]. Homer developed a model for medical technology adoption based on demandside (user dispositions to accept or abandon based on social exposure and evaluation of product performance) and supply-side (R&D for product performance improvement and investment in promotional activities) dynamics [33]. These models are extensions of the Bass Diffusion generic structure that includes a word-of-mouth diffusion process (social exposure and imitation) and external adoption from advertising, which enhance the realism of innovation diffusion [21].

In their systematic review of SD models on innovation systems, Uriona and Grobbelaar point to a "promising stream of research" based on Technological Innovation Systems (TIS) theory, that departs from the innovator-imitator structure of Bass Diffusion [34] (p. 34). TIS theory posits that the formation of a new technological innovation system requires seven key interacting elements: (1) entrepreneurial activities, (2) knowledge development, (3) knowledge diffusion, (4) guidance of search, (5) market formation, (6) mobilization of resources, and (7) creation of legitimacy [35–37]. The complex interactions of these elements determine the growth prospects of a new technology. Subsequently, Walrave and Raven operationalized the theory into a conceptual SD simulation model [36,38]. The main advantage of the TIS framework is its explanatory power for the market formation of new technologies—a "complex non-linear interactive process" that involves several actors and institutions [34] (p. 28). Indeed, market formation requires collective market-oriented action to develop "shared market infrastructure" for "supporting the functioning of a stable market" [39] (p. 244).

Our work contributes to existing knowledge by synthesizing the TIS model with medical technology adoption models. Similar to Paich, Peck and Valant, we represent the patient flow of the health care and prosthetics care system for amputees. Like Homer, we represent the demand-side dispositions of amputees to adopt or abandon the emerging digital prosthetic device based on performance evaluation and word-of-mouth diffusion. We then captured the supply-side conditions of digital prosthetics using the TIS framework. In doing so, we present a feedback perspective to digital prosthetics market formation prospects and the corollary effects on product adoption as well as on the socio-economic and health indicators for the amputee population.

#### *2.2. Data Collection*

The iterative model building process, from conceptualization, quantification, to validation, was conducted in collaboration with Toyota Mobility Foundation—expert in system dynamics and human centered design for promoting mobility (second author) and ProsFit Technologies—digital prosthetics service provider (third author). It is tradition in SD to include problem owners and experts in the model building process who possess important domain expertise, experiential knowledge, and mental models of the system under investigation [40–42]. During this process, several iterations of the model were presented to the collaborators for validation. In terms of model parameterization, ProsFit provided numerical estimates for some parameter values where existing data was not available. In such instances, ProsFit relied on its network of prosthetists and other experts in the field to corroborate their assumptions and understanding. Estimates and comments from these domain experts were anonymized and shared via email correspondence. Such estimates represent the best available data and expert judgement at the time of the model development.

Apart from expert opinion, existing peer-reviewed literature was utilized extensively for model conceptualization—especially so for the conceptual market formation subsystem in the model. As for quantification, parameter values were obtained either from epidemiological data reported in the literature or from secondary datasets (Table 1). The model is calibrated to data from the United Kingdom (UK), since expert knowledge and literature pertaining to the country is more readily available, but it can nevertheless be calibrated to other contexts.

**Table 1.** Data sources used for model parameterization.

#### *2.3. Model Description*

In this section, we present a simplified stock-and-flow representation of the model. The simplified structure is split into the top-level Health Care System (Figure 1) and the Market Subsystem (Figure 2). We then briefly describe the key feedback processes involved (see Appendix A for a more detailed description of each feedback loop).

*Systems* **2023**, *11*, x FOR PEER REVIEW 6 of 30

**Figure 1.** Simplified feedback structure of the amputee health care system. The sign on each link (+/–) indicates the polarity of the relationship between the connected two variables. **Figure 1.** Simplified feedback structure of the amputee health care system. The sign on each link (+/−) indicates the polarity of the relationship between the connected two variables.

Market Subsystem (Figure 2). We then briefly describe the key feedback processes in‐

volved (see Appendix A for a more detailed description of each feedback loop).

**Figure 2.** Simplified feedback structure of market subsystem. The sign on each link (+/–) indicates the polarity of the relationship between the two connected variables. **Figure 2.** Simplified feedback structure of market subsystem. The sign on each link (+/−) indicates the polarity of the relationship between the two connected variables.

The health care system for amputees is represented as an aging chain that captures the flow of people across different stages or compartments. This structure is further ar‐ rayed to better capture and represent the characteristics and choices of different popula‐

Amputees then flow into the prosthetic care stage, either into traditional prosthetics or digital prosthetics depending on the respective market share, from the primary care stage. They may achieve full mobility if successfully fitted with a prosthesis. However, both prosthesis types need to be replaced every three years on average [59]. As more am‐ putees enter the prosthetics care stage, prosthesis degradation over time increases the number of amputees awaiting replacement of their devices before re‐entering the pros‐ thetic fitting process again. In this regard, the *Prosthesis Lifecycle* loop (R1) could result in a growing pressure for prosthetics demand, emanating from our best efforts to success‐ fully fit new amputees with a prosthesis. Alternatively, amputees may dropout from the prosthesis fitting process altogether or abandon the device due to an unsuccessful fit

models, can help identify accumulations and key bottle‐necks in patient flows [27]. Figure 1 captures the flow of people from the general population and the PAD population stocks to acute care for either trauma‐related or PAD‐related amputation. At this stage, the *Pre‐ vention Pressure* loop (B1) works to reduce the PAD incidence rate over the long‐term as the PAD‐related amputations increase over time. This balancing feedback loop represents the prevention pressure faced by public health agencies to address the prevalence of PAD by stepping up efforts towards primary prevention, including early screening, smoking cessation, nutritional and activity programs [57,58]. A decline in PAD incidence would lead to a reduction of the PAD prevalence over time, which would eventually decrease

major lower‐limb amputations from PAD.

2.3.1. Health Care System

#### 2.3.1. Health Care System

The health care system for amputees is represented as an aging chain that captures the flow of people across different stages or compartments. This structure is further arrayed to better capture and represent the characteristics and choices of different population groups (e.g., age and prosthesis type). Aging chains, commonly used in SD health models, can help identify accumulations and key bottle-necks in patient flows [27]. Figure 1 captures the flow of people from the general population and the PAD population stocks to acute care for either trauma-related or PAD-related amputation. At this stage, the *Prevention Pressure* loop (B1) works to reduce the PAD incidence rate over the long-term as the PAD-related amputations increase over time. This balancing feedback loop represents the prevention pressure faced by public health agencies to address the prevalence of PAD by stepping up efforts towards primary prevention, including early screening, smoking cessation, nutritional and activity programs [57,58]. A decline in PAD incidence would lead to a reduction of the PAD prevalence over time, which would eventually decrease major lower-limb amputations from PAD.

Amputees then flow into the prosthetic care stage, either into traditional prosthetics or digital prosthetics depending on the respective market share, from the primary care stage. They may achieve full mobility if successfully fitted with a prosthesis. However, both prosthesis types need to be replaced every three years on average [59]. As more amputees enter the prosthetics care stage, prosthesis degradation over time increases the number of amputees awaiting replacement of their devices before re-entering the prosthetic fitting process again. In this regard, the *Prosthesis Lifecycle* loop (R1) could result in a growing pressure for prosthetics demand, emanating from our best efforts to successfully fit new amputees with a prosthesis. Alternatively, amputees may dropout from the prosthesis fitting process altogether or abandon the device due to an unsuccessful fit [16,60]. These individuals flow into the limited mobility stock. However, amputees with limited mobility may later decide to readopt a prosthesis and therefore re-enter the prosthetic care stage again. This process is captured in the *Prosthetics Re-entry* loop (R2). Both loops engender a reinforcing mechanism that moves amputees through different stages of prosthetics care. They do not independently multiply the number of amputees in the loop beyond those already within the closed aging chain.

With the introduction of a digital prosthetics market, amputees are probabilistically referred to a digital prosthetist dependent on the market share. The perceived success of digital prosthetics is then conceptualized as the ratio of the rate of successful digital fitting relative to traditional fitting. When the rate of digital fittings surpass the incumbent traditional technology, we can expect a stronger favorable word-of-mouth diffusion about the success or reputation of digital prosthesis [61,62]. Over time, we expect the reputation of digital technology to reinforce the growth of the digital market size and thus the market share of the digital prosthetics through the *Digital Growth* loops (R3 and R4). With a higher market share, even more amputees are more likely to be referred to a digital prosthetist or may seek out one themselves if they are re-adoptees. Concurrently, the *Prosthesis Attractiveness* loops (R8 and R9) encourage stronger uptake of digital prosthesis devices. Word-of-mouth diffusion about the relative success of digital technology could motivate individuals to stick to the process and thus translate to a lower drop-out rate. It could also motivate those who have previously abandoned the process to re-enter the fitting the process, consequently increasing the re-adoption rate. Here, the diffusion processes are driven by evaluations of the relative performance of digital prosthetics (successful digital fitting rates vs. traditional), similar to the adoption structure in Homer's model [33].

While the digital growth and prosthesis attractiveness loops drive the accumulation of amputees in the digital prosthetics care sector, the *Access Constraint* loops (B2, B3 and B4) counteract their reinforcing effects. Amputees' access to the prosthetics care stage is limited by the capacity of the sector (number of fittings that can be accommodated by available prosthetists). Fitting demand is driven by new amputees, those seeking to replace their degraded device, and re-adoptees who previously abandoned the fitting process. When

the fitting demand outweighs the fitting capacity, prosthetic accessibility reduces and thus limits amputees from entering the stage even if so desired. Over the longer term, however, the *Market Access* loops (R5, R6 and R7) work to improve accessibility for digital prosthetics by expanding capacity. A higher market share of digital technology would lead to an expansion of digital prosthetic clinics, which enables the sector to accommodate a larger number of patients. This effect, however, is delayed as it takes time to assess the market and set up new clinics.

As for the incumbent traditional prosthetics, they compete with the growing digital prosthetics market for referrals from the primary care stage. This is captured in the patient flow between the two stages, where entry to traditional prosthetics is determined by the inverse of the digital market share (i.e., 1–"Digital Market Share"). However, in the model, traditional prosthetics is assumed to be unaffected by the diffusion processes of digital technology. Should digital technology gain dominance, one might expect the incumbent's reputation to diminish and, as a result, more amputees might dropout and fewer might re-adopt a traditional prosthesis. Moreover, traditional prosthetics sector might also face a capacity contraction. Yet, these effects were not modeled for a more conservative estimate of digital prosthetics' impact since, without concrete data, this could add to further uncertainty of the model. Instead, the patient flows within the traditional prosthetics care sector were held at constant fractional rates, apart from new amputee referrals.

### 2.3.2. Market Subsystem

In the top-level health care system, we sought to explain the effects of digital prosthetics market growth on the prosthetics care sector. The complexity involved in market formation within the market subsystem is represented in Figure 2.

The process of technological knowledge development is described by the synergistic interaction of the *Technology Development* loop (R10) and *Knowledge Diffusion* loop (R11). Innovation development and diffusion of knowledge is required for any TIS to grow, and this is dependent on the level of resources available for R&D [35,63]. As innovation is developed and diffused through the exchange of knowledge between various actors in the system, the guidance of search for the new technology increases. Guidance of search refers to the "visibility and clarity" of the state of the art [35] (p. 423) that reflects the "promises and expectations of the emerging technology" [63] (p. 56). This helps in the priority-setting of that technology and directing more resources for further R&D, which would enable even more technological knowledge development and diffusion.

This process, in turn, attracts new entrepreneurs into the emerging market through the *Innovation Attractiveness* loop (R12) and *Knowledge Attractiveness* loop (R13). Entrepreneurs are central to any TIS for carrying out market-oriented action [63]. As more innovation is developed and diffused, the technological legitimacy of the technology increases and accumulates the perceived legitimacy of innovation system. [36]. This encourages more entrants to enter the market and grow the level of market-oriented entrepreneurial activity. Since entrepreneurial activities indicate the health and sustainability of an innovation system [35], this would bring in more external funding/resource stream into the system from private or public actors [36,63]. External funding further reinforces the growth of entrepreneurial activities through the *External Engine* loop (R14). External backing reduces the perceived entrepreneurial risks involved, and consequently is better able to attract further entry into the market [36,63]. Moreover, the external funding stream increases the total resources available in the system, which spurs more development of innovation that increases the legitimacy of the technology even further.

While the external engine stimulates entrepreneurial activity initially, the *System Legitimacy* loop (R17) endogenously generates internal ("financial, material, human capital") resources over the longer term for market sustainability [63] (p. 57). This loop comprises the two smaller *Internal Engine* loop (R15) and the *System Building* loop (R16), and is capable of driving the entire system [37]. Entrepreneurs contribute to the "development of formal market rules, establishment of intermediary networks, the building of infrastructure, or the

development of formal regulations" [36] (p. 1837). The developed market infrastructure generates market legitimacy for the TIS, which reduces market formation uncertainty and the perceived cost to participation [39]. Hence, more entrepreneurs are willing to overcome perceived risks and enter the market, contributing to further infrastructure development (R16). Moreover, these established market structures, mediated by entrepreneurial activities, "contribute to the creation of a demand for the emerging technology" [63] (p. 56). This increases the market size for the technology that generates internal resources from the market (R15). The synergy of the loops, reflected in R17, thus drives the self-reinforcing growth of entrepreneurial activities, market infrastructure, market legitimacy, and market size.

The formation of a niche market for the emerging technology, however, precipitates "resistance from actors with interests in the incumbent" regime [63] (p. 57). For instance, "when regime actors try to influence public discourses, or lobby against favourable support" [36] (p. 1837). This process is captured in the *Resistance* loop (R18). Regime resistance decreases the market legitimacy of the emerging technology, which disincentivizes entrants due to higher perceived risks. In turn, there will be less market infrastructure development to counter the regime resistance. As the niche market grows and competes with the incumbent regime, resistance could also come in the form of innovation. Given the new threat, regime actors would "increase their efforts to improve the performance of the existing regime through innovation" [36] (p. 1838). This is referred to as the sailing ship effect [64,65] and is represented in the *Sailing Ship* loop (B7). It contributes to a stronger regime resistance and counteracts the effects of the System Legitimacy loop.

Finally, the top-level health care system is connected to market subsystem through the *Digital Growth* loops and *Market Access* loops. As the reputation of digital fittings grow, we expect it to bolster the technological legitimacy of digital prosthetics. This would lend strength to the System Legitimacy loop, which ultimately increases the market size. With a larger market size, the market share of digital prosthetics rises relative to the incumbent traditional prosthetics. Importantly, the number of digital prosthetic clinics also increases to expand the fitting capacity. This improves the digital prosthetics accessibility, which enables more amputees to be fitted with a prosthesis and achieve mobility.

#### *2.4. Model Validation*

The described feedback structure was operationalized into a SD simulation model. The model was built in Stella Architect version 3.0 (SD modelling software from isee systems) using Euler Integration with a time-step of 1/16 of a month, or about 2 days, which is less than half of the smallest time constant of 7 days (0.23 months) for the pre-operation hospital stay in the primary care sector. The model is simulated over a time horizon of 480 months, representing January 2010 to January 2050. Simulation modelling facilitates the visualization of the impact of digital prosthetics on the health care system and, more importantly, experimentations to better understand the dynamic complexity of the system [21,29]. Here, we summarize the results of the model validation procedure as proposed by Forrester and Senge [66] and Barlas [67] to build confidence in the simulation results. A more detailed validation report is available in a previous iteration of this work [30].

The model structure is supported by relevant literature and input from stakeholders. As a digital prosthetics service provider (ProsFit) and a double lower-limb amputee, the third author of this paper was heavily consulted during the iterative process of model building to validate the structures in the health care system. Parameterization of the health care system was based on empirical data sources (Table 1). In instances where data was not available, the values were estimated from expert opinion. This pertains to the fractional dropout and readoption rates, which are estimates from ProsFit and their network of prosthetists. Parameter verification for the market subsystem, however, was challenging given the conceptual nature of the model. Thus, the parameter values set in the original model [38] was kept and subject to further sensitivity tests. All parameters and variables in the model were assigned units of measurement that are both mathematically

and conceptually consistent. The model documentation provided in the Supplementary Materials details the above for each variable in the model.

Moreover, direct- and indirect-extreme conditions tests were performed to ensure robustness of structural formulations. There were no computational errors detected in the model and the results conform to values that are within bounds. We further conducted sensitivity analysis for all parameters in the model. Each parameter was varied over 100 sensitivity runs. The variation is based on a uniform distribution random draw using Sobol Sequence sampling method [68]. The results of the sensitivity analysis are summarized in Table 2. Expectedly, the model was mostly sensitive to parameters in the conceptual market subsystem. As a result, it introduces uncertainty to the relatively empirical top-level model. This means that this model cannot produce accurate numerical estimations. Nevertheless, it was deemed more useful to represent the complexity of market formation than the alternative: a simplistic table function with high levels of sensitivity. Understanding how digital prosthetic market formation may plausibly occur from a feedback perspective could be useful to decision-makers seeking to improve mobility outcomes and maximize the impact of limited resources.

**Table 2.** Parameters resulting in model sensitivity.

\* Refer to Appendix B for the confidence plots of the model's sensitivity.

#### **3. Simulation Results**

#### *3.1. Baseline Setup*

The model was initialized in equilibrium to produce the baseline simulation results. In a previous iteration, we attempted to initialize the stocks at the obtained or calculated initial values [30]. However, there are virtually no numerical estimates for individuals in the various transitory stages of the primary care continuum and prosthetic care continuum. Consequently, we opted to initialize the stocks in their long-term equilibrium values to prevent transient stock adjustments. Moreover, initializing the model in equilibrium enables us to observe the full effects of any shocks exogenously introduced to the model—in our case, the formation of a niche digital prosthetics market.

To set the model in equilibrium, we held the total population of the UK constant at about 61.1 million individuals over the time horizon and initialized the population stocks in their long-term equilibrium values. The equilibrium switch in the market subsystems initializes the innovation diffusion stocks at zero and cuts off the exogenous input of relative resources. This ensures there are no dynamics in the market subsystem, thus representing a scenario wherein the prosthetics sector is solely serviced by the traditional prosthetics service providers at their existing capacity.

#### *3.2. Baseline Results*

The baseline results provide the estimated equilibrium values of the respective key indicators for the system (see Table 3). With a constant population size of 61.1 M people, we estimate a total of 84.8 K major lower-limb amputees, with about 85% of them being deemed medically eligible for a prosthetic device (71.9 K people). Of those eligible, only 5.5 K amputees are estimated to be fitted with a (traditional) prosthesis, thus resulting in a mobility proportion of over 7%. The amputee mobility proportion represents the proportion of eligible amputees who have achieved full mobility through a successful prosthesis fitting, which is determined by three factors. First, the accessibility of prosthetics services, which represents the percentage of demand that is met by the existing service capacity. The model estimates this to be just under 12%, indicating a bottleneck in the prosthetic care system. This is also within range of the WHO global estimate—that only 5 to 15% of amputees have access to prostheses [1]. Second, individual dispositions to drop out from the fitting process (estimated by experts to be about 10% for each stage of the process prior to the final device) and to readopt it (about 20%) at a later time. Third, the probability of final device fit success, which is about 50% for traditional devices [15,16].

**Table 3.** Baseline results of key indicators.

Furthermore, health-related economic indicators were calculated based on the inputs from the endogenous processes in the health and prosthetics care systems. The monthly economic productivity of amputees is estimated to be about USD14 M per month. This indicator represents the economic participation of amputees from returning to work and reintegrating into the workforce. It is conceptualized as the product of the estimated number of employed amputees with the gross domestic product per capita. Amputees not fitted with prostheses are excluded from workforce participation given their limited mobility. This is a gross simplification that does not completely reflect the economic contribution of amputees from other measures such as individual consumption. The total economic cost incurred, on the other hand, is estimated to be USD210 M per month. This includes the differentiated health care costs, unemployment and social payments, and the opportunity costs borne by families for caretaking. Lastly, the estimated total cost of prosthesis provision is about USD1.94 M per month, which includes both successful and failed prosthesis fittings. Reimbursements for prosthesis costs are assumed to be fully covered by national insurance mechanisms, which would otherwise be borne as out-of-pocket payments.

#### *3.3. Experimental Setup*

To simulate and investigate the impact of digital prosthetics on the baseline behavior, we introduced dynamics in the market subsystem module. This was done by setting the parameter value of Relative External Resources (RER) Size above 0 from month 96 for an assumed duration of 180 months (year 2018 to 2033) to exogenously kick start the dynamics. This simulates the deployment of external funding streams to support initial market growth

prior to self-sufficiency. Moreover, the Innovation Developed and Knowledge Diffused stocks were initialized at 0.01 to push them out of unstable equilibrium.

As mentioned, given the conceptual nature of the market subsystem, we have introduced several assumptions in the parameter values that would result in estimation errors. Therefore, we ran a global sensitivity analysis with combined variations in all the parameters that the model is sensitive to (as identified in Table 2). The experimental results, in turn, show the confidence intervals (up to 95%) of the key indicators from 1000 runs based on Sobol Sequence sampling method [68]. A total of 1000 runs was sufficient to fully explore the state space of the stocks in the market subsystems. This experiment gives us the full range of possibilities for digital prosthetics market growth and thus enables us to observe the corollary effects on the more empirical prosthetics care system.

## *3.4. Experimental Results*

Given its known sensitivity, the model produces a range of market growth for digital prosthetics, from 0.03% to 96% market share by 2050 with a mean of 43.6% (Figure 3). In general, with a RER size of more than 0, the External Engine loop (R14) powers the endogenous market formation processes that allows the Digital Market Share to start growing. However, as R14 loop is cut off by year 2033, we observe three behavioral patterns: (1) steady decline, (2) a much slower albeit continued growth, or (3) sustained growth. The sustainability of market growth is ultimately dependent on the strength of the Internal Engine (R15) and the System Legitimacy (R17) loops in endogenously generating sufficient internal market resources. *Systems* **2023**, *11*, x FOR PEER REVIEW 12 of 30 (1) steady decline, (2) a much slower albeit continued growth, or (3) sustained growth. The sustainability of market growth is ultimately dependent on the strength of the Inter‐ nal Engine (R15) and the System Legitimacy (R17) loops in endogenously generating suf‐ ficient internal market resources.

**Figure 3.** Confidence plot of the market share of the digital prosthetics service. **Figure 3.** Confidence plot of the market share of the digital prosthetics service.

Nevertheless, we can observe the impact of digital prosthetics market growth on the amputee population and mobility outcomes. As digital prostheses are introduced to the prosthetics care system, Figure 4a shows an increase in the mobility proportion from the baseline of 7% to a mean of 23% by 2050 (range: 6.2–50%). The introduction of digital prosthetics not only increases the existing service capacity, but also results in more suc‐ cessful fittings—a synergistic product of the Digital Growth (R3 and R4) and Prosthetics Attractiveness (R8 and R9) loops. The Digital Growth loops enable more amputees to en‐ ter the prosthetics care system either as a new entrant or a re‐adoptee, whereas the Pros‐ thetics Attractiveness loops discourages amputees in the fitting process from dropping out and encourages previous dropouts to re‐adopt a device. In turn, more amputees achieve mobility. The increased mobility further leads to improved health outcomes, in‐ cluding a lower mortality risk. In this sense, the expansion of digital prosthetics prevents more deaths, which accounts for the increase in the amputee population from the baseline of 84.8 K to an average of 87.6 K individuals by 2050 (range: 83.2 K–94 K) as seen in Figure Nevertheless, we can observe the impact of digital prosthetics market growth on the amputee population and mobility outcomes. As digital prostheses are introduced to the prosthetics care system, Figure 4a shows an increase in the mobility proportion from the baseline of 7% to a mean of 23% by 2050 (range: 6.2–50%). The introduction of digital prosthetics not only increases the existing service capacity, but also results in more successful fittings—a synergistic product of the Digital Growth (R3 and R4) and Prosthetics Attractiveness (R8 and R9) loops. The Digital Growth loops enable more amputees to enter the prosthetics care system either as a new entrant or a re-adoptee, whereas the Prosthetics Attractiveness loops discourages amputees in the fitting process from dropping out and encourages previous dropouts to re-adopt a device. In turn, more amputees achieve mobility. The increased mobility further leads to improved health outcomes, including a lower mortality risk. In this sense, the expansion of digital prosthetics prevents more deaths, which accounts for the increase in the amputee population from the baseline of 84.8 K to an average of 87.6 K individuals by 2050 (range: 83.2 K–94 K) as seen in Figure 4b.

Additionally, we observe that the mobility proportion develops similarly to the dig‐ ital market share; if the market share were to decline after year 2033, so would the mobility proportion. However, the mobility proportion develops at slower pace and changes to a

fitting process and (2) the effect of the Access Constraint loops (B2, B3 and B4) from pros‐ thetics accessibility that limits the number of amputees according to available service ca‐

4b.

pacity.

As for the prosthetics accessibility, it increases for the first of half the simulation du‐ ration before declining again, which consequently limits the growth of the mobility pro‐ portion. With reference to Figure 5, the peak of the mean accessibility is about 44% some time in 2028 (range: 13.6% to 100%), which eventually declines to 20% by 2050 (range: 2.8% to 93.6%). In general, the accessibility ratio increases as the digital market grows and adds additional service capacity to the existing level. In conditions where there is limited mar‐ ket growth, we observe that the accessibility declines around the time when the exogenous Additionally, we observe that the mobility proportion develops similarly to the digital market share; if the market share were to decline after year 2033, so would the mobility proportion. However, the mobility proportion develops at slower pace and changes to a smaller extent. This is due to (1) the multiple delays involved in the aging chains of the fitting process and (2) the effect of the Access Constraint loops (B2, B3 and B4) from prosthetics accessibility that limits the number of amputees according to available service capacity.

funding is cut‐off in year 2033. However, under more optimistic market growth condi‐ tions, we observe that the accessibility peaks prior to the cut‐off time and declines there‐ after. This is due to the higher volume of demand for replacing degraded prostheses gen‐ erated by the Prosthesis Lifecycle loop (R1). Prostheses have a lifecycle of 3 years on aver‐ age, and hence there is a captive consumer base that will continue to shore up fitting de‐ mand—more so when the proportion of fitted amputees is high. As seen in Figure 5, the total accessibility may increase to the maximum (100%) in instances where parameters enable a rapid and large expansion of digital clinics (e.g., Sensitivity of Clinics to Market Size) to meet the demand for fittings. Even then, it declines by the tail end of the simula‐ tion for the reasons described. As for the prosthetics accessibility, it increases for the first of half the simulation duration before declining again, which consequently limits the growth of the mobility proportion. With reference to Figure 5, the peak of the mean accessibility is about 44% some time in 2028 (range: 13.6% to 100%), which eventually declines to 20% by 2050 (range: 2.8% to 93.6%). In general, the accessibility ratio increases as the digital market grows and adds additional service capacity to the existing level. In conditions where there is limited market growth, we observe that the accessibility declines around the time when the exogenous funding is cut-off in year 2033. However, under more optimistic market growth conditions, we observe that the accessibility peaks prior to the cut-off time and declines thereafter. This is due to the higher volume of demand for replacing degraded prostheses generated by the Prosthesis Lifecycle loop (R1). Prostheses have a lifecycle of 3 years on average, and hence there is a captive consumer base that will continue to shore up fitting demand—more so when the proportion of fitted amputees is high. As seen in Figure 5, the total accessibility may increase to the maximum (100%) in instances where parameters enable a rapid and large expansion of digital clinics (e.g., Sensitivity of Clinics to Market Size) to meet the demand for fittings. Even then, it declines by the tail end of the simulation for the reasons described.

**Figure 5.** Confidence plot of the total accessibility of prosthetics services in terms of ability to meet total demand. As a result of the developments in the prosthetics care system, we can further assess As a result of the developments in the prosthetics care system, we can further assess the impact on the health-related economic indicators. The economic productivity of amputees follows a similar development to the mobility proportion since employed amputees make up a fraction of those who are mobile. Figure 6a shows and increase in the monthly productivity of amputees from an average of USD14 M to USD43 M by year 2050 (range: USD11.5 M to USD97.7 M). Again, these figures are underestimates that only partially captures the true economic contribution of amputees. Whereas Figure 6b shows a reduction in the monthly economic cost incurred, decreasing from USD210 M to USD202 M on average (range: USD289 M to USD211 M). The economic cost per capita reduces as the mobility proportion increases because mobile amputees incur smaller health care costs, social payments, and opportunity costs for their families. However, note that effect from the per capita cost reduction has been counteracted by the overall increase in amputee

the impact on the health‐related economic indicators. The economic productivity of

population size from improved health outcomes (less deaths). In this sense, the reduction in total economic cost is not as pronounced as the per capita reduction in economic cost. Based on these figures, we can further anticipate the net benefit of digital prosthetics service provision: the sum of the additional economic productivity and the amount of reduction in economic cost. The average net social benefit is then calculated to be a mean of USD37 M per month. mand—more so when the proportion of fitted amputees is high. As seen in Figure 5, the total accessibility may increase to the maximum (100%) in instances where parameters enable a rapid and large expansion of digital clinics (e.g., Sensitivity of Clinics to Market Size) to meet the demand for fittings. Even then, it declines by the tail end of the simula‐ tion for the reasons described. *Systems* **2023**, *11*, x FOR PEER REVIEW 14 of 30

**Figure 4.** (**a**) Confidence plot of the proportion of eligible amputees who are fitted with a prosthesis;

As for the prosthetics accessibility, it increases for the first of half the simulation du‐ ration before declining again, which consequently limits the growth of the mobility pro‐ portion. With reference to Figure 5, the peak of the mean accessibility is about 44% some time in 2028 (range: 13.6% to 100%), which eventually declines to 20% by 2050 (range: 2.8% to 93.6%). In general, the accessibility ratio increases as the digital market grows and adds additional service capacity to the existing level. In conditions where there is limited mar‐ ket growth, we observe that the accessibility declines around the time when the exogenous funding is cut‐off in year 2033. However, under more optimistic market growth condi‐ tions, we observe that the accessibility peaks prior to the cut‐off time and declines there‐ after. This is due to the higher volume of demand for replacing degraded prostheses gen‐ erated by the Prosthesis Lifecycle loop (R1). Prostheses have a lifecycle of 3 years on aver‐ age, and hence there is a captive consumer base that will continue to shore up fitting de‐

*Systems* **2023**, *11*, x FOR PEER REVIEW 13 of 30

(**a**) (**b**)

(**b**) Confidence plot of the total amputee population.

**Figure 5.** Confidence plot of the total accessibility of prosthetics services in terms of ability to meet total demand. **Figure 5.** Confidence plot of the total accessibility of prosthetics services in terms of ability to meet total demand. reduction in economic cost. The average net social benefit is then calculated to be a mean of USD37 M per month.

#### *3.5. Scenario Setup 3.5. Scenario Setup*

The experimental results have shown that there are three fundamental behavior modes within the range of possibilities for digital market growth. To better visualize the differentiated impact on the prosthetics care system, we developed three hypothetical sce‐ narios for digital market growth. First, a pessimistic scenario to represent the growth and steady decline. Second, a realistic scenario wherein the digital market experiences a much slower rate of growth after external funding is cut off. Third, an optimistic scenario to The experimental results have shown that there are three fundamental behavior modes within the range of possibilities for digital market growth. To better visualize the differentiated impact on the prosthetics care system, we developed three hypothetical scenarios for digital market growth. First, a pessimistic scenario to represent the growth and steady decline. Second, a realistic scenario wherein the digital market experiences a much slower rate of growth after external funding is cut off. Third, an optimistic scenario to represent the sustained market growth throughout the simulation duration.

represent the sustained market growth throughout the simulation duration. To this end, we conducted a sensitivity analysis with only the parameters that the model is behaviorally sensitive to (see Table 2) for a total of 50 runs. From these runs, we To this end, we conducted a sensitivity analysis with only the parameters that the model is behaviorally sensitive to (see Table 2) for a total of 50 runs. From these runs, we

selected the set of parameter values that produced the appropriate behavior mode for each

selected the set of parameter values that produced the appropriate behavior mode for each of the scenarios. These values are reported in Table 4.

**Table 4.** Behaviorally sensitive parameters and corresponding values for each scenario.

The set of parameter values for the pessimistic scenario results in a condition where there is a low level of resources flowing in the market subsystem. With a relatively lower RER size, the External Engine loop (R14) has a weaker reinforcing effect in the pessimistic scenario as compared to the other two. The weight of entrepreneurial activity modulates the level of market infrastructure development. A higher weight implies that infrastructure development is more dependent on the level of entrepreneurial activity in the system than the volume of resources available for market formation. Not only is there a low level of resources to begin with, but the market infrastructure development is also not as reactive to those resources in the pessimistic scenario. On the other end, in the optimistic scenario there are ample of resources in the system for market formation. A relatively lower weight on entrepreneurial activity further implies that market development is stimulated by the resources available. The realistic scenario represents a more likely median between the two extremes.

#### *3.6. Scenario Results*

We reproduced the three behavior patterns representing the varied conditions for digital prosthetics market growth (see Figure 7a). Under pessimistic market growth conditions, the market share of digital prosthetics growths to a maximum of about 5% before declining to 0.5% by 2050. Under the realistic scenario, the digital market share increases to 36% in 2033 and thereafter increases gradually to 43% by 2050. Whereas digital prosthetics experiences sustained growth in the optimistic scenario, capturing 80% of the market share by 2050.

Based on these three hypothetical market growth scenarios, we can observe the relative impact on the amputee mobility outcome in Figure 7b. In general, the mobility proportion follows the same behavioral pattern as the digital market share. The proportion increases as market share increases and vice versa. The gap between the realistic and optimistic scenarios for amputee mobility is disproportionately smaller than the gap for the digital prosthetics market share. This is due to the dampening effect of the Access Constraint loops (B2, B3 and B4) as explained before. Amputee mobility is being constrained by the fitting capacity that is unable to meet the demand. By further expanding digital fitting capacity, we can strengthen the effect of the Market Access loops (R5 and R6) to better counteract the constraint loops. In this sense, we can anticipate an even higher mobility proportion in the optimistic scenario than the 32% mobility observed in Figure 7b.

under three scenarios.

*Systems* **2023**, *11*, x FOR PEER REVIEW 16 of 30

**Figure 7.** (**a**) Comparative graphs of the digital prosthetics market share under three growth scenar‐ ios; (**b**) Comparative graphs of the proportion of eligible amputees who are fitted with a prosthesis under three scenarios. **Figure 7.** (**a**) Comparative graphs of the digital prosthetics market share under three growth scenarios; (**b**) Comparative graphs of the proportion of eligible amputees who are fitted with a prosthesis under three scenarios. (**a**) (**b**) **Figure 7.** (**a**) Comparative graphs of the digital prosthetics market share under three growth scenar‐ ios; (**b**) Comparative graphs of the proportion of eligible amputees who are fitted with a prosthesis

better counteract the constraint loops. In this sense, we can anticipate an even higher mo‐ bility proportion in the optimistic scenario than the 32% mobility observed in Figure 7b.

better counteract the constraint loops. In this sense, we can anticipate an even higher mo‐

As for the health‐related economic indicators, we can now graphically represent the net benefit of digital prosthetics service provision compared to the baseline behavior (Fig‐ ure 8a). In all three scenarios, the introduction of digital prosthetics results in a positive net benefit. A 0.5% digital market share in 2050 still yields a net benefit of USD5 M per month in the pessimistic scenario. This figure is USD41 M and USD53 M for the realistic and optimistic scenario, respectively. Moreover, we can compare the scale of the net ben‐ efit to that of the additional prosthesis reimbursement (Figure 8b). The additional reim‐ bursement is the difference between the total costs for prosthesis services and the baseline costs. In contexts where prosthetics services are covered by national health care systems, digital market growth increases the total costs borne by the state in terms of insurance reimbursements as the volume fittings increases over time. Nevertheless, additional costs incurred in that instance is far outweighed by the net benefit accrued—on average by a factor of 15 across all scenarios. As for the health-related economic indicators, we can now graphically represent the net benefit of digital prosthetics service provision compared to the baseline behavior (Figure 8a). In all three scenarios, the introduction of digital prosthetics results in a positive net benefit. A 0.5% digital market share in 2050 still yields a net benefit of USD5 M per month in the pessimistic scenario. This figure is USD41 M and USD53 M for the realistic and optimistic scenario, respectively. Moreover, we can compare the scale of the net benefit to that of the additional prosthesis reimbursement (Figure 8b). The additional reimbursement is the difference between the total costs for prosthesis services and the baseline costs. In contexts where prosthetics services are covered by national health care systems, digital market growth increases the total costs borne by the state in terms of insurance reimbursements as the volume fittings increases over time. Nevertheless, additional costs incurred in that instance is far outweighed by the net benefit accrued—on average by a factor of 15 across all scenarios. As for the health‐related economic indicators, we can now graphically represent the net benefit of digital prosthetics service provision compared to the baseline behavior (Fig‐ ure 8a). In all three scenarios, the introduction of digital prosthetics results in a positive net benefit. A 0.5% digital market share in 2050 still yields a net benefit of USD5 M per month in the pessimistic scenario. This figure is USD41 M and USD53 M for the realistic and optimistic scenario, respectively. Moreover, we can compare the scale of the net ben‐ efit to that of the additional prosthesis reimbursement (Figure 8b). The additional reim‐ bursement is the difference between the total costs for prosthesis services and the baseline costs. In contexts where prosthetics services are covered by national health care systems, digital market growth increases the total costs borne by the state in terms of insurance reimbursements as the volume fittings increases over time. Nevertheless, additional costs incurred in that instance is far outweighed by the net benefit accrued—on average by a factor of 15 across all scenarios.

ment cost incurred under three scenarios. **Figure 8.** (**a**) Comparative graphs of the undiscounted net benefit of digital prosthetics under three growth scenarios; (**b**) Comparative graphs of the undiscounted additional prosthesis reimburse‐ ment cost incurred under three scenarios. **Figure 8.** (**a**) Comparative graphs of the undiscounted net benefit of digital prosthetics under three growth scenarios; (**b**) Comparative graphs of the undiscounted additional prosthesis reimbursement cost incurred under three scenarios.

#### **4. Discussion**

In summary, our simulation model allows exploring the range of mobility outcomes for the amputee population given different market growth conditions for digital prosthetics. We observed in all experimental scenarios that an increase in digital prosthetics market share was associated with improved mobility outcomes. Specifically, there was an increase

in the mobility proportion, an increase in total economic productivity, a decrease in total economic costs and an overall positive net benefit from digital prosthetics even under pessimistic conditions. Our model could, thus, serve as a tool for health policy planners to explore a shift in prosthetics service provision to better mitigate the accessibility problem and bolster the social impact of prosthesis usage.