Practical Use of Regularization in Individualizing a Mathematical Model of Cardiovascular Hemodynamics Using Scarce Data

Individualizing physiological models to a patient can enable patient-specific monitoring and treatment in critical care environments. However, this task often presents a unique “practical identifiability” challenge due to the conflict between model complexity and data scarcity. Regularization provides an established framework to cope with this conflict by compensating for data scarcity with prior knowledge. However, regularization has not been widely pursued in individualizing physiological models to facilitate patient-specific critical care. Thus, the goal of this work is to garner potentially generalizable insight into the practical use of regularization in individualizing a complex physiological model using scarce data by investigating its effect in a clinically significant critical care case study of blood volume kinetics and cardiovascular hemodynamics in hemorrhage and circulatory resuscitation. We construct a population-average model as prior knowledge and individualize the physiological model via regularization to illustrate that regularization can be effective in individualizing a physiological model to learn salient individual-specific characteristics (resulting in the goodness of fit to individual-specific data) while restricting unnecessary deviations from the population-average model (achieving practical identifiability). We also illustrate that regularization yields parsimonious individualization of only sensitive parameters as well as adequate physiological plausibility and relevance in predicting internal physiological states.


INTRODUCTION
Clinical care automation has been a domain of interest for a few decades by virtue of its potential for error-free and vigilant performance of routine and low-level patient monitoring and treatment tasks (Hemmerling et al., 2010;Salinas et al., 2011;Dussaussoy et al., 2014;Rinehart et al., 2015;Brogi et al., 2017;Hundeshagen et al., 2017;Pasin et al., 2017), yet realizing this potential is contingent upon establishing the safety and the performance characteristics of clinical care automation. Patient physiology models built upon physical principles (hereafter called physiological models) can facilitate the development (Jin-Oh et al., 2012;Bighamian et al., 2014;Jin et al., 2018) and the testing (Kovatchev et al., 2009;Ortiz et al., 2010;Brown et al., 2015) of clinical care automation capabilities. However, individualizing physiological models (which can enable the systematic development and the testing of patient-specific clinical care automation) presents formidable challenge due to the conflict between model complexity and data scarcity. That is, physiological models are complex and involve a large number of unknown patient-specific parameters that frequently exhibit interaction properties, whereas clinical data pertaining to an individual patient are scarce in both quantity and quality. For example, blood pressure (BP) can increase in response to an increase in cardiac output (CO) and an increase in total peripheral resistance (TPR), but routine clinical measurements often provide only BP which does not offer sufficient information to uniquely determine CO-and TPR-related model parameters.
In fact, such a conflict appears to be a main obstacle in characterizing a wide range of mechanistic models in physiology as well as in broader fields of biology and physics, which has been called the sloppiness property (Transtrum et al., 2010(Transtrum et al., , 2015Machta et al., 2013;White et al., 2016) or lack of practical identifiability (Raue et al., 2009;Maiwald et al., 2011) in different contexts. If not properly addressed, such conflict may yield a physiological model suffering from non-unique and physically irrelevant parameter values as well as poor individual-specific internal (i.e., unmeasured) state predictions.
The existing body of work has addressed the practical identifiability challenge mainly in two ways: (i) by modifying the model structure or parameterization toward more simplified or lumped models and (ii) by utilizing some form of regularization (or more generally, prior knowledge) to additionally inform the modeling procedure and effectively reduce the complexity of the model search space. Analyzing and modifying the model structure has particularly strong precedents in physiological modeling, where effective tools such as variation-based sensitivity analysis (Sobol, 2001;Saltelli et al., 2008Saltelli et al., , 2010 and profilelikelihood (Raue et al., 2009(Raue et al., , 2014Maiwald et al., 2016) methods are used to assess the appropriateness of a model structure and the importance of model parameters with respect to given data. While useful in informing potential modifications to the model structure and parameterization (Machta et al., 2013;Mattingly et al., 2018), many of these tools are formalized on a predetermined level of data availability, making them less powerful if data availability varies (which is quite common in clinical scenarios). Regularization is an established framework to cope with conflicts between model complexity and data scarcity, by which one can potentially compensate for the scarce data available to individualize a physiological model based on prior knowledge. To date, the notion of prior distribution (Toni et al., 2009;Gelman et al., 2013) and its closely related counterpart in regularization (Benning and Burger, 2018) have been pursued in a wide range of deterministic and stochastic modeling problems to incorporate additional information into the modeling problem and control the complexity of the model search space (Davidian and Giltinan, 2003;Tivay et al., 2019;Zhang et al., 2019).
The goal of this work is to garner potentially generalizable insight into the role of regularization in individualizing a complex physiological model using scarce data by studying its effect in a clinically significant critical care case study of blood volume (BV) kinetics and cardiovascular hemodynamics in hemorrhage and circulatory resuscitation. We construct a population-average physiological model as prior knowledge and individualize the physiological model via regularization to elucidate how regularization enables the physiological model to capture salient individual-specific characteristics (leading to the goodness of fit to individual-specific data) while restricting unnecessary deviations from prior knowledge (achieving practical identifiability). We also analyze how regularization affects the plausibility and the relevance of the physiological model under varying levels of scarcity in individualspecific data in terms of its parameter values and internal state predictions.
This paper is organized as follows: section "Cardiovascular Hemodynamics in Hemorrhage and Circulatory Resuscitation" outlines the physiological model and the experimental data used in the case study. Section "Materials and Methods" describes a step-by-step procedure for applying regularization to individualize the physiological model. Section "Results" presents the results, which are discussed in section "Discussion". Section "Conclusions" concludes the paper with suggested future work.

CARDIOVASCULAR HEMODYNAMICS IN HEMORRHAGE AND CIRCULATORY RESUSCITATION
Credible physiological models of BV kinetics and cardiovascular hemodynamics have the potential to facilitate the development and the testing of circulatory resuscitation devices and algorithms . Prior work has reported a wide spectrum of physiological models (in terms of scale and detail) to represent volume kinetics (Hahn, 2010;Bighamian et al., 2016) and cardiovascular hemodynamics (Moss et al., 2012;Bighamian et al., 2018) in response to stimuli. In some cases, these physiological models were also analyzed for practical identifiability (Moss et al., 2012;Marquis et al., 2018;Pathmanathan et al., 2019;Pironet et al., 2019). In this work, we employ an extended and enhanced version of a physiological model reported in our prior work , which concerns the effect of hemorrhage and circulatory resuscitation on the changes in BV kinetics and cardiovascular hemodynamics.

Model Description
The physiological model is composed of component mathematical models that represent BV kinetics, cardiovascular function, and BP regulation, all at the systems level (Figure 1). The BV kinetics model is a lumped-parameter model that simulates the effect of fluid gain, blood loss, and intercompartmental fluid exchange on changes in BV. In this model, the time rate of change in BV ( V B ) can be described by the following differential equation: where u, v, and u ex , respectively, denote the rate of external fluid gain (e.g., fluid resuscitation), external fluid loss (e.g., hemorrhage and urinary output), and internal fluid exchange between BV and interstitial tissue compartments. The rate of fluid exchange is macroscopically modeled as a control input signal that regulates BV to a target value r BV as follows: where K p is the controller gain and r BV represents the target change in BV, which is computed according to the history of external fluid gain and loss as follows: where the ratio 1 1+a u determines the fraction of the total external fluid gain that eventually remains within the BV compartment (with the rest shifted to the interstitial tissue compartment), and the ratio 1 1+a v determines the fraction of the total fluid loss that is not compensated for by a refill from the interstitial tissue compartment. In our prior work, we demonstrated that the macroscopic expressions above can reproduce BV kinetic responses to volume perturbations with crystalloid, colloid, and blood (Bighamian et al., 2016. The cardiovascular function model is built to compute CO from BV based on the CO-venous return equilibrium principle (Guyton et al., 1975) and the left ventricular (LV) pressurevolume relationship (Lankhaar et al., 2009). The details of such a derivation can be found in our previous work , which yields the following model relating BV to CO: where HR (t) is the heart rate, TPR(t) is the TPR associated with the arterial circulation, C s is the systemic capacitance, E s is the left ventricular elastance, A and B define the shape of the LV pressure-volume loop, γ represents the approximate proportionality constant between central venous (P CV ) and LV end diastolic (P LVED ) pressures (γ ≈ P LVED /P CV ), k denotes the ratio between the resistance to venous return and TPR (k = R VR /TPR), V B0 is the initial BV, and η represents a proportional relationship between V B0 and the unstressed BV V BU : The BP regulation component of the model approximates the regulatory mechanism in the body that monitors the discrepancy between a pre-specified, intrinsic mean arterial BP (MAP) set point (MAP target ) versus the subject's actual MAP (MAP(t)) and adjusts TPR (and accordingly, the resistance to venous return) to maintain MAP near the set point level. This component is modeled with the following equations: where TPR 0 is the initial value of TPR, which can be found from the initial value of CO and MAP (TPR 0 = MAP 0 /CO 0 ), k TPR is the gain of MAP regulation, and p TPR represents the speed of MAP regulation. Due to the relatively short time period (180 min) associated with the dataset as well as the presumably maximal level of short-term BP regulation responses due to large initial hemorrhage, the full mechanistic details of BP regulation mechanisms (e.g., the long-term aspects associated with the responses of the renin-angiotensin-aldosterone and anti-diuretic hormone systems) are not considered, and the collective shortterm effects of BP regulation mechanisms are macroscopically expressed by Eq. 5.

Experimental Dataset
The experimental measurements acquired from animals in an array of prior work (Rafie et al., 2004;Vaid et al., 2006;Marques et al., 2017) were used as data in this case study. A total of N = 23 sheep constituted the population considered in this study. The data included the time series sequences of BV, CO, and MAP acquired from these sheep subjects. In each sheep, initial BV was measured using a dye concentration method. Then, subsequent changes in BV were measured with hemoglobin dilution after correction for the loss of red blood cells. CO was measured with a flow probe placed on the pulmonary artery. MAP was measured with arterial catheterization. These signals were measured at 5min intervals during a measurement period of 180 min. Each sheep was subjected to a large initial hemorrhage (25 ml/kg; see "H" at 0-15 min in Figure 4) followed by resuscitation (see "I" at 30-180 min in Figure 4) and subsequent minor hemorrhages (5 ml/kg; see "H" at 50-55 and 70-75 min in Figure 4). Fluid resuscitation was performed by previously developed closed-loop control algorithms (Rafie et al., 2004;Vaid et al., 2006;Marques et al., 2017) designed to maintain MAP. Hemorrhage and fluid resuscitation yielded diverse BV, CO, and MAP responses in the datasets, which is appropriate to analyze the role of regularization in individualizing the physiological model ( Table 1).

MATERIALS AND METHODS
In this section, we present a procedure for applying regularization to the physiological model outlined in section "Cardiovascular Hemodynamics in Hemorrhage and Circulatory Resuscitation". Special emphasis is given to the construction of prior knowledge, which involves the choice of a regularization function employed in individualizing the physiological model. Then, data analysis details are presented.

Individualizing the Physiological Model
Given the data associated with an individual subject i, the inputs given to the subject (i.e., hemorrhage and fluid resuscitation) are denoted by U d i = {u (t) , v(t)}, while the outputs of the subject are denoted by a vector Y d i . For instance, in case BV, CO, and MAP data are available in the individual i, then Y d i is defined as: where X d i t 1 , . . . , t n X is a row vector containing the measurements of X at times t 1 , . . . , t n X . The nominal noise covariance matrix associated with Y d To obtain model predictions Y i (θ, U d i ) corresponding to the elements of the data vector Y d i , the model equations in section "Cardiovascular Hemodynamics in Hemorrhage and Circulatory Resuscitation" are numerically solved given the inputs U d i and the vector of lumped model parameters θ, which is defined as follows (see Appendix Table A1 for details and ranges): (7) The problem of individualizing the physiological model can be formulated as solving the following maximum a posterioi estimation problem (Gelman et al., 2013): whereθ i is the vector of model parameters associated with the individual i, J 1 (θ) = − log P Y d i |U d i , θ corresponds to the likelihood of the individual-specific data with respect to the model output, and J 2 (θ) = − log P(θ) is a regularization term that corresponds to prior knowledge P (θ) about the parameter values. Assuming that P Y d i |U d i , θ is a Gaussian distribution due to noise on the output measurements, we can set J 1 (θ) as the covariance-weighted mean-squared output errors: In case there is no informative prior knowledge about the probability distribution of the model parameters, J 2 (θ) in Eq. 8 does not depend on θ and is effectively removed. However, in that case, the optimization problem in Eq. 8 is likely to become ill-posed and suffer from practical identifiability challenges when the physiological model is too complex relative to the available individual-specific data (which is frequently the case).
The main idea to address the practical identifiability problem in the absence of sufficient individual-specific data is to construct an appropriate prior distribution P(θ) to additionally inform the modeling problem. This prior distribution can be constructed based on the prior knowledge about plausible parameter ranges as well as heterogeneous data from a population of individuals. In this work, we formulate a prior in the shape of a Laplace distribution of the following form: where θ m is the m-th element in θ,θ m is the m-th element of the mode of the distributionθ, λ is the regularization weight, and b m is a normalization factor for θ m , which reflects prior knowledge about reasonable and/or physiological ranges for the model parameter [i.e., given physiological upper θ m and lower θ m bounds on θ m , we set b m = (θ m −θ m ); see Appendix Table A1 Frontiers in Physiology | www.frontiersin.org forθ m andθ m associated with the physiological model in section "Cardiovascular Hemodynamics in Hemorrhage and Circulatory Resuscitation"]. Equation 10 is a form of L 1 -regularization Jenatton et al., 2011), which represents the inter-individual variability observed in the population data by compressed parametric deviations from a mode denoted byθ. This representation may be possible if the studied physiological modeling problem has the sloppiness property (Machta et al., 2013;Transtrum et al., 2015), i.e., a wide range of parameter sensitivities, implying compressibility (see section "Role of Regularization in Individualizing Physiological Models" for more details as far as our case study is concerned). Having Eq. 10, the problem of constructing the prior reduces to determining the parameters λ andθ from heterogeneous population data. Techniques exist in the field of non-linear mixed-effects modeling (Davidian and Giltinan, 2003) and system identification (Pan et al., 2018) that can deal with such problems, often with comparable results (Kataria et al., 1994;Hahn et al., 2011). In this work, we defineθ as the solution to the following maximum-likelihood optimization problem based on the population data from L = N − 1 individuals: (11) Note that the population-based modeling problem in Eq. 11 is supposedly less susceptible to practical identifiability challenges than the individualization problem in Eq. 8 by virtue of the larger amount of heterogeneous data available from the population. The resulting model withθ will hereafter be called the populationaverage model, which is intended to represent an aggregated and sensible model of typical physiological behavior.
The regularization weight λ serves as the relative scaling factor between the spread of the prior in Eq. 10 and the spread of the likelihood in Eq. 9. To find an appropriate individual-specific weight λ * i , we adopt an approach related to the L-curve method (Hansen, 1992;Hansen and O'Leary, 1993). The weight λ * i may be estimated by solving the optimization problem in Eq. 8 for a range of λ ∈ [0 λ max ] and plotting the resulting values for the goodness of fit cost J 1 (θ) and the deviation distance D(θ) = J 2 (θ)/λ (Figure 2). An initial increase of λ from zero tends to decrease D by compressing parametric deviations fromθ in low-sensitivity directions. This increases J 1 , but the extent is not large relative to the decrease in D since the effect of restricting low-sensitivity (i.e., sloppy) deviations on J 1 is, by definition, small. This trend continues until λ reaches a critical value. However, beyond this critical value, sensitive deviations also become restricted, leading to a substantial increase in J 1 relative to the decrease in D. The appropriate λ * i can be chosen as the value of λ in between these regime changes (Figure 2).
In sum, individualizing a physiological model via regularization can be performed in two steps: (i) constructing a prior (regularization function) from heterogeneous data measured in a population by solving Eq. 11 and (ii) individualizing the physiological model with regularization using scarce individual-specific data by solving Eq. 8 (Figure 3).

Data Analysis
We applied the individualization procedure outlined in section "Individualizing the Physiological Model" to the physiological model outlined in section "Cardiovascular Hemodynamics in Hemorrhage and Circulatory Resuscitation". In particular, we constructed an appropriate regularization function and individualized the physiological model under varying levels of scarcity in the individual-specific data. To enable the data analysis with relatively small sample size, we performed a leave-one-out type of data analysis (details follow).
We individualized the physiological model under three levels of increasing data scarcity: (i) case 1, in which the physiological model was individualized given BV, CO, and MAP data, (ii) case 2, in which the physiological model was individualized given CO and MAP data, and (iii) case 3, in which the model was individualized given only MAP data. We considered these data scarcity scenarios for two investigational reasons: (i) comparison of the physiological models individualized with and without the use of regularization in the three cases may offer insight into how regularization affects the parameter values in response to varying data scarcity and (ii) analysis of the prediction errors for unmeasured internal states (i.e., BV in case 2 and BV and CO in case 3) with and without the use of regularization may offer insight into how regularization affects the physiological plausibility and the relevance of the physiological model under data scarcity.
For each animal subject, we analyzed the data pertaining to the rest of the L = 22 animal subjects with the maximumlikelihood method in Eq. 11 to derive a population-average physiological model, characterized byθ, that predicts BV, CO, and MAP responses to hemorrhage and fluid resuscitation inputs. Then, we individualized the physiological model to the (excluded) animal via regularization by solving Eq. 8 with Eq. 10 as the regularization function. We selected the FIGURE 3 | Individualizing physiological models with scarce data via regularization. In the first step, a prior (regularization) is constructed from heterogeneous data measured in a population by solving Eq. 11. In the second step, the physiological model is individualized with regularization using scarce data by solving Eq. 8. appropriate value of λ for each individual subject as described in section "Individualizing the Physiological Model" (Figure 2). We repeated individualizing the physiological model in this way for the three dataset scarcity cases. For comparative analysis, we also individualized the physiological model to each animal without using regularization for all the three data scarcity scenarios.
With the goal of assessing the plausibility and the relevance of the physiological models, we compared the models individualized with and without regularization as well as the populationaverage model. We devised two quantitative error metrics for this purpose: (i) output prediction error e 1 , defined as the normalized root-mean-square error (RMSE) between model-predicted and measured output signals (among BV, CO, and MAP) explicitly used for individualizing the physiological model and (ii) state prediction error e 2 , defined as the normalized RMSE between model-predicted and measured internal state signals not used in individualizing the physiological model. For example, in case 3, e 1 is the normalized RMSE associated with MAP: while e 2 is the normalized RMSE associated with BV and CO: We employed the Wilcoxon signed-rank test in conjunction with the Bonferroni correction to assess the statistical significance of the difference in these metrics between the models. A numerical simulation of the physiological model was performed using MATLAB R 's ODE solvers in the Simulink R environment. The numerical optimization was performed using MATLAB R 's Optimization Toolbox. Data analysis and visualization was performed using the seaborn and the matplotlib libraries in Python. Table 2 shows the output (e 1 ) and the internal state (e 2 ) prediction errors associated with the population-average model and the individualized models (both with and without the use of regularization), while Figure 4 presents the measured versus the model-predicted output and internal state signals associated with the population-average and the individualized physiological models. Figure 5 presents the representative behaviors of the output prediction error (e 1 ), the internal state prediction error (e 2 ), and the parametric deviation distance (D) with respect to the regularization weight λ in an individual subject. Figure 6 presents the regularized parametric deviations from the population-average model with respect to varying degrees of data scarcity.

DISCUSSION
With the long-term goal of advancing the future development and the testing of clinical care automation based on physiological models, we sought to investigate the role of regularization in individualizing a physiological model using scarce data.
Our findings from the case study of blood volume kinetics and cardiovascular hemodynamics in hemorrhage and circulatory resuscitation suggest that regularization can be effective in individualizing the physiological model to capture salient individual-specific characteristics and restrict unnecessary deviations from the population-average model. Furthermore, regularization results in appropriately varying parametric deviations to cope with the varying levels of data scarcity, thereby securing the physiological plausibility and the relevance of the individualized physiological model (details follow).

Role of Regularization in Individualizing Physiological Models
First, the use of regularization effectively maintained the goodness of fit to individual-specific data while reducing parametric deviations in insensitive directions. In all scenarios associated with the three levels of data scarcity (see section "Data Analysis"), the difference in e 1 was small between the physiological models individualized with and without the use of regularization (Table 2 and Figure 4). In contrast, the parametric deviations [in terms of the distance D(θ) = J 2 (θ)/λ averaged across all the animals] were smaller when regularization was employed (1.19 in case 1, 1.14 in case 2, and 0.30 in case 3) than when it was not (1.87 in case 1, 1.89 in case 2, and 2.82 in case 3). This suggests that regularization assists in individualizing a physiological model to achieve an adequate level of goodness of fit while restricting parametric deviations that may only result in overfitting.
Second, the physiological model individualized with the use of regularization exhibited improved physiological plausibility and relevance in comparison with the one individualized without the use of regularization in terms of the accuracy in predicting internal state signals. In general, the goodness of fit associated with BV and CO responses was deteriorated when the corresponding measurement was removed in individualizing the physiological model. Figure 4 shows that indeed (i) BV prediction is deteriorated when BV data become unavailable (case 1 → case 2 and case 3) and (ii) CO prediction is likewise deteriorated when CO data become unavailable (case 2 → case 3). In case 2, e 2 (i.e., BV prediction error) was smaller when regularization was employed than when it was not employed (p = 0.03) ( Table 2). In case 3, e 2 (i.e., BV and CO prediction errors) was likewise smaller when regularization was employed than when it was not employed (p < 0.01) ( Table 2). Interestingly, MAP prediction was notably improved as data became more scarce when regularization was not employed [case 1 (7.23 mmHg) to case 2 (6.94 mmHg) to case 3 (4.85 mmHg), all in terms of mean errors; see Table 2]. Thus, the large deterioration in e 2 may be attributed to overfitting to MAP data. In conjunction with the finding above on the role of regularization in restricting unnecessary parametric deviations, this suggests that regularization improves the ability of the individualized physiological model to predict the internal states accurately by preventing the undesired over-tuning of model parameters away from their population-average values. One additional minor note is that BV prediction achieved with regularization was comparable between case 2 and case 3. This can be interpreted as follows: (i) both case 2 and case 3 did not use BV data in individualizing the physiological model, and regularization tended to reduce the BV component in the physiological model to the population-average model (as can be seen by the comparable BV prediction errors associated with these models; see Table 2) and (ii) presumably CO data did not provide much implications on the behavior of BV on top of the prior knowledge [as can be seen by the BV prediction errors in case 2 with and without the use of regularization (0.25 L and 0.29 L, both in terms of mean errors; see Table 2]. Third, the individualized physiological model exhibited goodness of fit parametric deviation behavior with respect to the regularization weight λ as anticipated in Figure 2, especially in the case of very scarce data (case 3; see Figure 5). A sharp initial decrease in D(θ) = J 2 (θ)/λ relative to a modest increase in e 1 was observed when λ was initially increased from zero. Further increase in λ after a critical point caused e 1 to largely increase relative to a steady decrease in the parametric deviation distance and eventually approach that of the population-average model. In sum, an appropriate λ * i could be selected as outlined in section "Individualizing the Physiological Model" (λ * = 0.15).
Interestingly, there appeared to be a notable link between the behaviors of D(θ) and e 2 : e 2 initially showed a notably decreasing trend similar to that of D(θ), which suggests that regularization improves the ability of the individualized model to predict internal states (Figure 5). For larger values of λ, e 2 increased, which can be attributed to the restriction of sensitive parametric deviations, degenerating the individualized model to the population-average model. Fourth, the physiological model individualized with the use of regularization was associated with a significantly superior output prediction error and a comparable internal state prediction error in comparison with the population-average model ( Table 2). This suggests that regularization is effective in judiciously individualizing sensitive model parameters to achieve desirable FIGURE 5 | Representative trends of the output prediction error e 1 , the internal state prediction error e 2 , and the parametric deviation distance in a subject with respect to the regularization weight λ.  Table A1). Vertical axis: model parameters.
goodness of fit to individual-specific data while at the same time capitalizing on prior knowledge (in the form of regularization) to preserve the ability to adequately predict unmeasured internal states.

Parametric Deviations in Relation to Varying Data Scarcity
In our case study, regularization-induced deviations from the population-average model showed two noteworthy trends. First, with regularization, deviations tended to decrease as the data became more scarce (Figure 6). In other words, (i) regularization allowed meaningful deviations if ample data were available so that the physiological model could be better individualized by absorbing the individual-specific characteristics represented by the data, whereas (ii) it restricted unwarranted deviations if scarce data were available so that the physiological model could fall back to the population-average model. This can be viewed as a desirable behavior in that increasingly leveraging the prior knowledge (i.e., the population-average model) as the data scarcity increases is the intended effect of regularization. In contrast, such a parametric deviation trend was not observed when regularization was not used. The opposite occurred instead: scarce data resulted in more aggressive deviations from the population-average model (see the deviation distance results reported in section "Role of Regularization in Individualizing Physiological Models"), which consequently compromised the fidelity of the individualized physiological models (e.g., its internal state prediction was deteriorated; see Table 2). Second, regularization appeared to effectively individualize those parameters relevant to the presented data. For example, the most visibly deviated parameter in case 3 was the one associated with BP regulation (k TPR ), which is reasonable in that only MAP measurements were presented to individualize the physiological model in case 3 ( Figure 6C). The parameters associated with both CV function (γ/BC s and η) and BP regulation (k TPR and MAP target ) were likewise visibly deviated in case 2, in which CO and MAP measurements were presented ( Figure 6B). Finally, the parameters associated with BV were likewise largely deviated (α v , in particular) when BV measurements were presented in addition to CO and MAP data ( Figure 6A). One caveat is that, despite the presence of these trends, the studied parameters belong to a complex physiological system and it may not be possible to attribute each parametric deviation completely to the availability of a specific type of measurement.

Study Limitations
This work has limitations. First, this work investigated the role of regularization in a specific case study. Hence, the insights obtained from the case study may not be universally true in all physiological modeling problems. In particular, the legitimacy of the key assumption used to address the adverse effect of scarce data on the quality of individualized models (that the inter-individual variability observed in the population data can be represented by compressed deviations from a population-average model according to the Laplace distribution in Eq. 10) may depend on the specificity of the physiological modeling problem at hand. On the one hand, such a compressed representation may indeed be valid in many real-world physiological modeling problems with the sloppiness property (Machta et al., 2013;Transtrum et al., 2015) (which implies parameter-space compressibility). On the other hand, regularization may not prove effective if the quality of the data and the model structure are such that the majority of the model parameters are associated with sensitive deviations (where individualizing the physiological model requires deviations in all the model parameters; in this case, regularization is not needed) or insensitive deviations (where a population-average model may suffice because of negligible inter-individual variability). In any case, the validity of the "parameter-space compressibility" assumption in a specific physiological modeling application can be assessed based on model prediction errors and parameter values in the individualized physiological models: the bias introduced by the inadequacy of the assumption will manifest itself as poor goodness of fit of the individualized physiological model to the data, while a lack of compressibility will manifest itself as estimated elements ofθ i in Eq. 8 that all deviate fromθ in Eq. 11. Second, regularization does not guarantee, and may even prevent, the convergence of the physiological model to the "true" individual-specific physiological model. Indeed its main purpose is to minimize the deviations from the prior knowledge (i.e., population-average model). Thus, the quality of a physiological model individualized with regularization, especially using highly scarce data, hinges upon the quality of prior knowledge. Hence, the physiological model individualized with regularization is a candidate approximation to the ideal individual-specific physiological model (which would be obtained with sufficiently diverse and high-quality data obtained from an individual). Taking these aspects into account, regularization may be viewed as an effective tool to individualize a physiological model in cases with a relatively complex model and scarce individualspecific data.

CONCLUSION
We sought to garner potentially generalizable insight into the role of regularization as prior knowledge in individualizing physiological models using scarce data. The findings obtained from a clinically significant case study illustrated that regularized individualization creates physiological models equipped with several desirable properties: (i) an adequate trade-off between goodness of fit to scarce individual-specific data and deviation from a population-average model, (ii) improved physiological plausibility and relevance, and (iii) parametric deviations relevant to the scarcity level of the data. Noting that regularization is not prevalently used in the physiological modeling domain, future work must be devoted to exploring the use of regularization in individualizing a range of physiological models and developing appropriate theory and novel structures for the regularization functions based on appropriate physiological assumptions.

DATA AVAILABILITY STATEMENT
The datasets analyzed in this study can be obtained from Dr. George C. Kramer at the University of Texas Medical Branch (gkramer@utmb.edu).

AUTHOR CONTRIBUTIONS
AT, CS, and J-OH formulated the procedure to individualize the physiological model via regularization and wrote and revised the manuscript. XJ, AL, CS, and J-OH developed the physiological model. AT and J-OH analyzed the data. XJ and AL reviewed the manuscript.

ACKNOWLEDGMENTS
We thank Dr. George C. Kramer at the University of Texas Medical Branch for provision of the experimental dataset for analysis as well as his physiological consultation and guidance in data analysis.