Sepsis Reconsidered: Identifying novel metrics for behavioral landscape characterization with a high-performance computing implementation of an agent-based model

Objectives: Sepsis affects nearly 1 million people in the United States per year, has a mortality rate of 28–50m% and requires more than $20 billion a year in hospital costs. Over a quarter century of research has not yielded a single reliable diagnostic test or a directed therapeutic agent for sepsis. Central to this insufficiency is the fact that sepsis remains a clinical/physiological diagnosis representing a multitude of molecularly heterogeneous pathological trajectories. Advances in computational capabilities offered by High Performance Computing (HPC) platforms call for an evolution in the investigation of sepsis to attempt to define the boundaries of traditional research (bench, clinical and computational) through the use of computational proxy models. We present a novel investigatory and analytical approach, derived from how HPC resources and simulation are used in the physical sciences, to identify the epistemic boundary conditions of the study of clinical sepsis via the use of a proxy agent-based model of systemic inflammation. Design: Current predictive models for sepsis use correlative methods are limited by patient heterogeneity and data sparseness. We address this issue by using an HPC version of a system-level validated agent-based model of sepsis, the Innate Immune Response ABM (IIRBM), as a proxy system in order to identify boundary conditions for the possible behavioral space for sepsis. We then apply advanced analysis derived from the study of Random Dynamical Systems (RDS) to identify novel means for characterizing system behavior and providing insight into the tractability of traditional investigatory methods. Results: The behavior space of the IIRABM was examined by simulating over 70 million sepsis patients for up to 90 days for the following parameters: cardio-respiratory-metabolic resilience; microbial invasiveness; microbial toxigenesis; and degree of nosocomial exposure. In addition to using established methods for describing parameter space, we developed two novel methods for characterizing the behavior of a RDS: Probabilistic Basins of Attraction (PBoA) and Stochastic Trajectory Analysis (STA). Computationally generated behavioral landscapes demonstrated attractor structures around stochastic regions of behavior that could be described in a complementary fashion through use of PBoA and STA. The stochasticity of the boundaries of the attractors highlights the challenge for correlative attempts to characterize and classify clinical sepsis. Conclusions: HPC simulations of models like the IIRABM can be used to generate approximations of the behavior space of sepsis to both establish “boundaries of futility” with respect to existing investigatory approaches and apply system engineering principles to investigate the general dynamic properties of sepsis to provide a pathway for developing control strategies. The issues that bedevil the study and treatment of sepsis, namely clinical data sparseness and inadequate experimental sampling of system behavior space, are fundamental to nearly all biomedical research, manifesting in the “Crisis of Reproducibility” at all levels. HPC-augmented simulation-based research offers an investigatory strategy more consistent with that seen in the physical sciences (which combine experiment, theory and simulation), and an opportunity to utilize the leading advances in HPC, namely deep machine learning and evolutionary computing, to form the basis of an iterative scientific process to meet the full promise of Precision Medicine (right drug, right patient, right time).

There are nearly 1 million cases of sepsis in the United States each year, with a mortality rate between 2 28-50% [1] While operational care process improvements in the last 20 years that have led to reduction in 3 mortality [2] therapeutic options for sepsis remain variations of anti-microbial and physiological support 4 dating back nearly a quarter century, and, crucially, there remains no FDA approved biologically 5 targeted therapeutic for the treatment of sepsis. In an era where the overall goal of medical care is to 6 provide personalized/precision medicine, which should mean "right drug for the right patient at the right 7 time," achieving this goal in sepsis is significantly hampered by the lack of success in translating basic 8 science knowledge into robust and effective therapeutics -a problem pervasive across the biomedical 9 spectrum [3]. Over a quarter century of failed attempts to modulate sepsis raises a key question: how 10 tractable are current investigatory strategies, namely increasing granular clinical data collection and 11 analysis coupled with pre-clinical experimental investigations, aimed at achieving this goal? Is there a 12 fundamental limitation of the capabilities of the current approaches to attempting to characterize and 13 control sepsis? This question of the epistemic limits of traditional research methodology is not limited to 14 sepsis, but is arising across the biomedical research spectrum. For example, a recent paper examines 15 such methodological limits by applying current analytical approaches used in neuroscience to a known, 16 simpler proxy information processing system, a microprocessor, and finds those approaches insufficient 17 for meaningful characterization of the behavior of the microprocessor [4]. We propose that a similar 18 approach can be used to evaluate the tractability of increasing reductionist attempts to define and 19 characterize clinical sepsis. We use a known computational model, a previously developed agent-based 20 model of sepsis recognized to be vastly simpler than the real world system, and apply to it the types of 21 data extraction and system state characterization possible with respect to clinical biomarker analysis. 22 This use of computational/simulation models as proxy systems has a long history in systems engineering, 23 and has been described specifically in the use of biomedical agent-based models (ABMs) [5]. The 24 rationale for using a simulation proxy model to establish analytical boundary conditions is based on the 25 fact that there is complete knowledge of the computational model, i.e. there are no "hidden variables," 26 and therefore it has an internally consistent ground truth. As a computer program, every component 27 and interaction of the model is known, including probabilities associated with stochastic events inputted 28 into the model's code. Given that the ABM is vastly simpler than the real world system, we assert that 29 if complete knowledge of the behavioral output of the simpler proxy model is insufficient to meaningfully 30 predict its behavior, attempts to do so in the real world will be similarly futile. We propose that an 31 iterative process utilizing large-scale simulation of sufficiently complex proxy models and advanced 32 analysis of simulation data can establish "boundaries of futility" with respect to what can or cannot be 33 known about clinical sepsis utilizing data-centric methods. This approach is only made possible by the 34 exponential growth of computing power seen in leadership-class high performance computing (HPC) that 35 now makes tractable the ability to near comprehensively characterize the behavioral landscape of an 36 ABM. We utilize an HPC implementation of a previously developed system/population-level validated 37 ABM of the innate immune response [6] (the Innate Immune Response ABM or IIRABM) to generate a 38 first approximation behavioral landscape of sepsis to guide the development of suitable metrics with 39 which to characterize that behavioral landscape. Though developed over 15 years ago, the IIRABM 40 remains valid insomuch it reproduces the overall dynamics of the acute inflammatory response, 41 producing recognizable system level outcomes (healing, pro-inflammatory death, hypo-immune death and 42 overwhelming initial infection) using a knowledge-based rule system that, while admittedly incomplete, 43 has not been demonstrated to be incorrect in its represented behaviors by accumulated discoveries since 44 its development. In fact, some of the behaviors generated by the IIRABM presaged their general 45 recognition within the sepsis community, specifically the temporal concurrence of pro-and 46 anti-inflammatory cytokine responses (as opposed to sequential pro-and compensatory responses) [7,8] 47 and the importance of the immuno-incompetent/attenuated recovery phase of sepsis, particularly with 48 respect to its prolonged duration [9][10][11]. For the current investigations, we acknowledge the 49 incompleteness of the IIRABM, considering it an abstract, substantially less complex representation of 50 the human innate immune response. As such, we further pose that the multi-dimensional behavioral 51 landscape of the IIRABM is less complex than can be expected of the real world system, and that 52 2/19 characterization of the behavioral landscape of the IIRABM represents a lower bound, first 53 approximation of the ability to characterize the real world system. The investigations presented here are 54 divided into two distinct but related tasks: 1) identification of the region of parameter space of the 55 IIRABM that corresponds to biologically/clinically plausible behavior, and 2) the development of novel 56 metrics and dynamical systems analysis to characterize the behavior of the IIRABM and provide insight 57 into the tractability and viability of current pathways toward precision medicine for sepsis.

59
For our investigations we use a previously validated agent-based model (ABM) of sepsis, the Innate 60 Immune Response ABM (IIRABM) [6]. The IIRABM was ported from NetLogo [12] to C++ and 61 implemented using MPI (Message Passing Interface) 2.0 for the parallelization and distribution of the 62 model. The IIRABM utilizes cells as its agent level. Individuals are represented as an interaction surface 63 between endothelial cells and circulating inflammatory cells and are subjected to variable perturbations 64 representing either infectious insult or tissue trauma. Multiple cell types are represented, including 65 endothelial cells, macrophages, neutrophils, TH0, TH1, and TH2 cells, as well as precursor cells for these 66 cell lines. Details for the IIRABM can be seen in Ref [6]. System damage is represented by an aggregate 67 measure of individual endothelial cell damage; the death threshold is set at 80%, representing the ability 68 of ICU care to render short-term lung, kidney, and cardiac failure to be survivable [6]. The increasing 69 death threshold is important for two reasons: it simulates ICU care in that modern medical techniques 70 allow a human to survive an amount of damage that would not have been survivable before modern 71 medicine; and an increasing death threshold allows the model to evolve to states which would have been 72 unreachable if execution were halted earlier [6,13]. Beyond the components and interactions of the 73 inflammatory response, the IIRABM utilizes 4 "external variables" that describe clinical features that 74 affect, but are not systemically integral to, the clinically-relevant systemic inflammation: environmental 75 toxicity (recurrent microbial exposure seen in a clinical environment), host resilience (an aggregate 76 variable representing the cardiorespiratory reserve of the patient and manifest as the ability of damaged 77 tissue to recover its oxygen level), and two measures of microbial virulence, invasiveness (ability to 78 spread in host tissue) and toxigenesis (ability to kill host tissue) (See Supplemental Table 1). Note that 79 these variables represent factors that clearly influence an individual's response during sepsis, but that 80 there are currently no clinically accessible, or even identifiable metrics that would allow categorization of 81 a particular patient. It is because of this combination of known influence and inability to categorize that 82 the entire parameter space needs to be characterized. In total, 8800 environmental parameter sets were 83 represented. A set of infectious injuries, the consequences of which ranged from trivial to severe, was 84 applied to each of these parameter sets. This injury was represented by a circular region with a radius 85 that varied from 1 to 40 grid-spaces in diameter. The inherent stochasticity built into the IIRABM 86 represents both the real/biological and epistemic variability seen in the immune response to injury and 87 infection; this allows the model to capture the variability seen in a clinical population. Thus, identical 88 initial injuries can lead to diverging patient trajectories depending on how the immune system of an in 89 silico patient responds to the simulated pathogen and injury distribution. The primary basis for the 90 behavioral validity of the IIRABM is its ability to reproduce the general dynamics of a system's 91 inflammatory response to infectious insult resulting in four clinically relevant trajectories: healing, 92 hyper-inflammatory system death, immuno-paralyzed late death, and overwhelming infection [6]. These 93 validity criteria are conserved from the initial description of the IIRABM, and are based on the common 94 sense observation that a biologically plausible system should not be always killed by the smallest insult 95 possible, nor should it always completely recover from the largest insult possible. Therefore, a 96 biologically plausible parameter set for the model would have a lower threshold of system injury from 97 which it would always recover (i.e. "healed") as well as an upper threshold of system injury from which 98 it would always die (i.e. "overwhelming infection). The additional dynamic classes of behavior become 99 evident within this zone of plausible behavior, and arise due to the fundamental purposes of the 100 inflammatory response, namely the eradication of infection by invasive microbes (enhanced by forward 101 feedback processes to increase its efficacy) and facilitation of the recovery from system damage (which 102 3/19 mandates negative feedback control to attenuate the initial response). Both of these dynamics result in 103 the effective eradication of invading microbes, but result in subsequent system death due to the 104 disordered internal processes. While each of these dynamic classes represents a different type of control 105 failure (pro-inflammatory death = distorted forward feedback and immune-paralysis = distorted negative 106 feedback) for purposes of this current initial behavioral analysis we classify them together as representing 107 internal system failure, clinically analogous to multi-organ failure. Future analysis will separate these 108 two dynamics, but within this current paper we define three possible system level outcomes: 1) healing, 109 2) system failure and 3) overwhelming infection. We use these three criteria as the basis of population technically feasible to represent this information in a single form. Therefore, we have developed as set of 119 nested metrics to emphasize specific perspectives aimed at establishing a baseline behavioral landscape 120 of clinical sepsis. We utilize five primary metrics with which to characterize the behavior of the system: 121 1) population-level outcome distributions for each parameter set 2) multi-dimensional parameter space 122 characterization across the entire range of parameter sets, 3) severity of outcome distribution heatmaps 123 to reify across biologically plausible parameter sets 4) probabilistic basins of attraction (PBoA) and 5) 124 stochastic trajectory analysis. The first 3 metrics are used to identify the region of parameter space of 125 the IIRABM that corresponds to biologically/clinically plausible behavior and confirm the baseline   Table 1). Each combination of these parameters was evaluated with N = 100 stochastic replicates for 131 each initial injury represented as a circular injury with a radius ranging from 1 to 40 cell widths. The results of the simulations (N = 100) for each parameter set shown as a histogram with injury size as 134 the x-axis, and the count of a specific outcome for the y-axis and is used to examine outcome-severity 135 progression for an individual parameter set as a function of the magnitude of the injury simulated. The 136 three possible outcomes of the IIRABM -healing, system failure (sepsis), or overwhelming infection, are 137 represented using three possible bars for each injury number. This approach was initially described in 138 rudimentary form in [6] and forms the basis by which the biologically/clinically plausible overall 139 parameter space characterization is performed. For the metrics listed below, Metric #2 can be 140 considered to utilize the collapse of Metric #1 into a 0-dimensional point data representation, whereas 141 Metric #3 can be considered to utilize the collapse of Metric #1 into a 1-dimensional line representation. 142 Figure 1 demonstrates 3 representative population outcome distributions as a function of injury size.

143
Where Figure 1A depicts a non-plausible condition where the system always dies with infection and 144 Figure 1B depicts a non-plausible condition where the system always recovers from infection, Figure 1C 145 depicts a plausible, clinically relevant condition where at some level the system heals, at some level the 146 system dies from infection, and in between there is the generation of behavior where the system clears Outcome counts for a specific injury size are shown on the y-axis. Blue shading indicates complete healing, green shading indicates death by sepsis, and yellow shading indicates death by infection. Panel A displays a histogram of outcomes as a function of injury size for simulations run with high bacterial virulence (invasiveness=4, toxigenesis=10), low host resilience (resilience=0.05), and high environmental toxicity (environmental toxicity=10); this parameter set is implausible because the system always dies with even the minimal injury. Panel B displays a histogram of outcomes as a function of injury size for simulations run with low bacterial virulence (invasiveness=1, toxigenesis=1), high host resilience (resilience=0.9), and low environmental toxicity (toxicity=1); this parameter set is implausible because the system always heals despite the maximal injury. In Panel C, the full range of outcomes as a function of initial injury size is displayed for an infection with moderate toxigenesis (toxigenesis=4) and invasiveness (invasiveness=2), low environmental toxicity (environmental toxicity=1), and average host resilience (resilience=0.1). Smaller initial injuries (shown on the left of the figure) more often result in complete healing of the simulated host (as expected) while the largest initial injury values, which also deliver the highest infectious load, result in death from overwhelming infection. The central region of this panel shows a region in which the host is capable of fighting off the applied infection, but is unable to recover from the inflammatory state that was required to fight the infection.

5/19
50% died due to system failure. We consider the population outcome distribution seen in 1C as an 151 example of a clinically relevant parameter set, since all 3 classes of outcomes (healing, system failure and 152 overwhelming infection) are present.

154
A multi-dimensional parameter sweep across a range of initial perturbation (initial infection) identifies 155 the upper and lower bounds with respect to clinically plausible behavior, as reflected in Metric #1 above. 156 The lower bound is that defined by parameter combinations where the system cannot be killed by the 157 minimal, least virulent microbes (i.e. always heals); the upper bound is defined by parameter 158 combinations where the system cannot recover from the maximal, most virulent microbes (always 159 overwhelmed by infection). Therefore, biologically plausible/clinically relevant parameter space resides 160 between these two boundaries. The representation of this parameter space is a multi-dimensional grid, 161 where each point in the grid represents a classification of the outcomes with respect to Metric #1; points 162 in that space that reflect biologically plausible behavior is then termed Clinically Relevant space. The Relevant space: while all 4 of the external variables represent "real" factors determining the severity of 193 sepsis, there is no current way of even qualitatively correlating those rankings with any clinically 194 accessible metric. This means that there is no way to know what the overall distribution is of all sepsis 195 patients across the Clinically Relevant space, and that any attempt to "match" a mortality rate from a 196 particular parameter set to that seen in clinical sepsis (i.e. the 20-50% mortality generally 197 reported [1,2]) would be a false and misleading attempt at demonstrating "validation."  simulations of a single parameter set. The shading indicates the percentage of deaths; blue represents no 218 deaths at and yellow represents 100% death. Figure 3A and Figure 3B show identical parameter sets but 219 different ordering, where patients represented in the Figure 3B received antibiotics twice a day for 10 220 days beginning 6 hours after their injury. As expected, the application of antibiotics greatly increases the 221 survivability of the simulated injuries -the yellow region is moved to the right, indicating improved 222 8/19 survival to higher initial injury numbers when antibiotics are available. The antibiotic parameter set was 223 reordered (again from worst to best) as the application of antibiotics changed the system dynamics. The 224 "stepping" feature that can be seen in Figure 3B is due to an emergent clustering of parameter sets by 225 the host-resilience parameter. The most resilient hosts are on the top and the least resilient are on the 226 bottom of each step group. This makes sense because the application of antibiotics significantly 227 decreases the bacterial infection's influence on the system; the other external parameters, invasiveness, 228 toxigenesis, and environmental toxicity, are all directly related to the simulated microbial infection.

229
When antibiotics are applied, the host-side parameters take precedence for outcome determination, a 230 plausible behavior further substantiating the validity of the IIRABM. Following the use of these first 3 231 metrics to identify the biologically plausible Clinically Relevant regions of parameter space for the model, 232 we then turn to the development of 2 novel data analysis methods in order to evaluate basic properties of 233 the system and how they might impact the investigation of real-world, clinical sepsis. The current approach to disease classification/prognosis/diagnosis is to use some set of measurements 236 (be they biomarkers or physiological signals) taken from the patient at a particular point in time, with 237 the goal that such information, or time series sequences of such information, will be able to distinguish 238 between groups with some threshold of acceptable success. The more general description of this concept 239 is characterizing the relationship between system state and system trajectory, the basis of the study of 240 Dynamical Systems [14] -that is, a system that evolves in time according to some rules in which future 241 states evolve from the current state. In dynamical systems theory, an attractor is a state or set of states 242 to which certain states tend to evolve [15]. A basin of attraction is the region of system state space zone without a clear structure across a range of system damage suggests that a particular level of either 274 of these cytokines alone poorly correlates with system outcome and would be an inadequate biomarker 275 for outcome (this is consistent with the current understanding of the utility of IL-10 or TNF as a marker 276 of disease severity). Figure 4C shows the probability of death as a joint function of TNF and IL10 levels. 277 Interestingly, the "slope" of the region of greatest mortality is not a constant, but levels off as the 278 concentration of IL10 increases. This suggests that an intervention augmenting systemic IL10 production 279 would be insufficient to significantly alter a patient's outcome or disease trajectory. Figure 4D, 280 demonstrating the PBoA for population levels of TH1 cells versus system damage, shows a qualitatively 281 similar structure as to Figure 4A and B: there appears to be virtually no independent predictive 282 capability for the number of TH1 cells present. Figure 4E displays the PBoA for population levels of TH2 283 cells against system damage, where the "slope" of the region of highest mortality has a slight upward 284 slope, suggesting that higher levels of TH2 cells can be associated with negative outcomes in sepsis.

285
However, the width of the vertical distribution of outcomes is large enough to wash out the predictive 286 capability of this trend. Alternatively, plotting the population of TH1 cells versus the population of TH2 287 cells ( Figure 4F) generates a more interesting PBoA that demonstrates the benefits of being able to shift 288 between multi-dimensional data projections. Here we see a very poorly bounded stochastic zone for lower 289 levels of both T-cell subtypes, likely representing the earlier phases of the disease course as the system 290 evolves, that transitions into a different region of state space with a more defined structure. In this 291 region the predictive capability of the TH1/TH2 ratio becomes more powerful, with an increase in the 292 proportion of TH2 cells associated with adverse outcomes, a finding that is consistent with both existing 293 reports [17] and the negative consequences of later immunoparalysis [9][10][11]. However, in general, the biomarkers for sepsis [18,19]. These PCA plots (Supplemental Material 1) show no discrimination at any 300 of the time points, consistent with the findings of the PBoAs. Given the predictive insufficiency of 301 system state snapshots, we now turn to more directly evaluating the trajectory space of the IIRABM. Reseeding Point Figure 6. This figure displays patient trajectories for a single parameter set (invasiveness=2, toxigene-sis=4, environmental toxicity=0, and host resilience=0.1) and single initial random number seed. The random number generator was re-seeded 100 times at at oxygen deficit=3000 (top), 3500 (middle) and 4000 (bottom). The original trajectory is shown in red and the trajectories generated from reseeding the random number generator are shown in blue.This image further reinforces the identification of a "stochastic zone" (a region of parameter space in which stochastic noise is (or can be) more powerful than the influence of either attractor).
has crossed out of the stochastic zone and into the basin of attraction for the Death Attractor.

353
The use of proxy models is ubiquitous and essential to the practice of science. Every laboratory model 354 (in vitro or in vivo) represents a proxy system for the real system being investigated. An essential aspect 355 of all proxy models is that they represent some degree of abstraction or divergence of detail with respect 356 to the real world system; in fact the breadth of their explanatory power is directly related to their 357 generalizability. However, for highly-engineered biological experimental models abstractions/limitations 358 are impossible to characterize formally, with a consequent impact on the ability to extrapolate findings 359 from those models to their reference systems and manifesting in the persistence of the Translational opposed to fitting a specific set of experimental/empirical data, we propose utilizing the ability of 366 simulation to vastly increase the "data coverage" of possible states and trajectories of the target system, 367 allowing a more comprehensive representation of the behavior of that system [20,21].

368
In Ref [22] the authors undertake the near monumental task of compiling the list of randomized 369 clinical trials that form the basis for current critical care practice (see Online Supplement to [22]). signaling theorem [23,24], which set lower bounds for a sufficient sampling criterion to discretize a 379 continuous signal. For example, if one were to try to reconstruct the signal y = sin(t) by sampling every 380 π seconds starting at t = 0, the reconstructed signal would simply be y = 0; all of the interesting features 381 of the function would be washed out by under-sampling. While the Nyquist-Shannon theorem deals with 382 1 dimensional functions we can extrapolate to the current sepsis work, and biomedical research in 383 general. In fact, the "Crisis of Reproducibility" noted in the biomedical arena [25] can be at least 384 partially attributed to the perpetual under-sampling of the total space of possible outcomes, with the 385 consequence that a misleading statistical significance for one particular study will become evident as 386 subsequent studies add to the portions of behavior space sampled [26].

387
If the sampling is perpetually temporally sparse (as will be the case for the foreseeable future due to 388 the logistical challenges of obtaining good quality clinical data), multiplication of dimensions across 389 which each time point is evaluated only increases the dimensions in which information becomes lost. As 390 we cannot derive a theoretically optimal sampling criterion (the model has no explicit functional form), 391 we instead perform a sufficiently large parameter sweep to bound the region of interest in parameter 392 space -in this specific case that region was the Clinically Relevant region in which any outcome was 393 possible -and provide a reference point by which the degree of under-sampling can be quantified.

394
Without some concept of the overall system behavior space that can approximate the scope of the 395 denominator problem in induction, the danger moving forward is that, given a specific data set, "some" 396 best fit function can be found for that specific data set, but will intrinsically be limited in it's 397 applicability to the more general condition. This explains why, in terms of clinical decision support as a 398 path towards precision medicine, physiology-based clinical screening tools [27][28][29] perform essentially as 399 well as multiplexed biomarker/-omics assays [30][31][32][33] for their respective predictive targets (onset of sepsis 400 and sepsis outcome). We assert that the perpetual under-sampling of sepsis behavior space places an 401

14/19
upper bound on the predictive capacity of any algorithm based on the under-representative data set 402 (which is a pragmatically fixed constraint). The methodology employed for the development of all these 403 scoring systems all assume that the data sets utilized sufficiently represent the overall set of behaviors 404 possible, which the current study demonstrates is not the case. In fact, given the "bowtie structure" of 405 biological signaling and regulatory networks [34], dynamic physiology-based prediction 406 systems [29] [35] [36] likely have a greater chance of success in refining individual patient trajectories, but 407 at the cost of not providing insight into the mechanistic drivers that would be potential targets for 408 therapeutic control. Future work with simulation-based characterization of sepsis will likely involve 409 mapping between computationally-generated behavioral landscapes and finer-grained temporal clinical 410 phenotype characterization using advanced physiological metrics to guide discovery of mechanistic 411 determinants of individual patient trajectories. The importance of identifying general properties of 412 patient trajectories, particularly with respect to the concept of attractors influencing those trajectories, 413 has a long history in the study of sepsis [37,38]. However, moving this recognition beyond the conceptual 414 level has been a challenge over the ensuing decades, due in main part to the fact that the concept of 415 "attractors" derive from dynamical systems theory, which is based on the assumption of deterministic 416 functions (even the study of stochastic dynamical systems primarily involve the addition of a noise term 417 to an existing deterministic function). It is readily clear that biological systems are not deterministic 418 (due to both real and epistemic stochasticity), and are more accurately described as Random Dynamical 419 Systems (RDS) [16]. As noted above, this recognition requires the modification of the classical attractor 420 concept to our Probabilistic Basins of Attraction (PBoA). The concept of the PBoA is predicated on the 421 fact that a "master equation" for the dynamical system is not known; the PBoA can only be derived 422 empirically. And since real-world biological data is too sparse, this requires the use of large-scale 423 simulation experiments to generate the landscapes of the PBoAs, which in turn is only possible with the 424 employment of HPC platforms. In time, for a sufficiently validated model, one can imagine the PBoA as 425 a valuable clinical diagnostic tool. Chemokine levels and their time-derivatives would be used as input 426 coordinates to a high-dimensional matrix containing outcome probabilities. These probabilities would 427 then inform future treatment paths. Since they are derived from mechanistic simulations, PBoAs can 428 also be used to generate targets and hypotheses for therapy discovery, e.g., if it were shown that the 429 death outcome is strongly associated with a specific signal configuration, this would suggest control interventions/control strategies guided by PBoAs and STAs to first determine controllability of the 439 overall system, and subsequently attempt to determine the robustness/scalability of control 440 strategies/policies identified. We propose that proxy simulation models, such as the IIRABM, can be 441 refined with an iterative loop between simulation and empiric plausibility evaluation (just as 442 meteorological simulations are constantly being refined), and during this process can be used to help 443 establish boundaries for plausible and implausible investigatory strategies, define expectations for 444 possible success and potentially eliminate futile approaches [41]. Absent this approach, currently, 445 prospective analysis of putative stratification systems (to say nothing of potential therapeutic 446 interventions) cannot be done from a pure real-world empirical, data-centric fashion: such failures will 447 only become evident ex post facto, with the consequent loss of time, money, resources and lives [42][43][44][45][46][47][48][49][50]. 448 Akin to the current use of HPC resources in physics and meteorology, the approach demonstrated in this 449 paper, the use of extremely large scale simulation and simulated data, provides the only demonstrated 450 effective scientific strategy to prospectively identify the boundaries of fruitful investigation. The use of 451 large-scale simulation based science for sepsis specifically, and biomedicine in general, represents a 452 potential path towards the full-scale application of engineering and control principles to the care of 453 15/19 individuals in terms of personalized/precision medicine and truly rational design of effective therapeutics. 454