{'fraction_non_alphanumeric': 0.04025753747100973, 'fraction_numerical': 0.01521340302537298, 'mean_word_length': 5.205455912361723, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 9, 'https://': 15, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mathematical Modeling, In-Human Evaluation and Analysis of Volume Kinetics and Kidney Function after Burn Injury and Resuscitation

Existing burn resuscitation protocols exhibit large variability in treatment efficacy. Hence, they must be further optimized based on comprehensive knowledge of burn pathophysiology. A physics-based mathematical model that can replicate physiological responses in diverse burn patients can serve as an attractive basis to perform non-clinical testing of burn resuscitation protocols and to expand knowledge on burn pathophysiology. We intend to develop, optimize, validate, and analyze a mathematical model to replicate physiological responses in burn patients. Using clinical datasets collected from 233 burn patients receiving burn resuscitation, we developed and validated a mathematical model applicable to computer-aided in-human burn resuscitation trial and knowledge expansion. Using the validated mathematical model, we examined possible physiological mechanisms responsible for the cohort-dependent differences in burn pathophysiology between younger versus older patients, female versus male patients, and patients with versus without inhalational injury. We demonstrated that the mathematical model could replicate physiological responses in burn patients associated with wide demographic characteristics and injury severity and that an increased inflammatory response to injury may be a key contributing factor in increasing the mortality risk of older patients and patients with inhalation injury via an increase in the fluid retention. The mathematical model may provide an attractive platform to conduct non-clinical testing of burn resuscitation protocols and test new hypotheses on burn pathophysiology.


INTRODUCTION
BURN is a leading cause of unintentional injury and death in the United States. According to a recent fact sheet from the American Burn Association, burn injury is the 8 th leading cause of death in adults >65 years old, 3 rd in children 5-9 years old, and 5 th in children 1-4 years old [1]. Each year, there are nearly 500,000 burn injury incidences only across the United States, among which 40,000 are hospitalized. Burn injury results in severe inflammatory responses that lead to increased leakage of intravascular water into the tissues, which, if not treated, engenders fatal consequences such as hypovolemic shock, ischemia, multiple organ failure, and generalized edema [2], [3]. The majority of burn patients treated in burn centers survive, but many of them suffer from complications [4], [5].
Currently, treatment of burn patients involves a resuscitation protocol in which the dose of the resuscitation fluid is titrated frequently to maintain the urinary output (UO) response in a burn patient to a therapeutic target range (e.g., 30-50 ml/hr or 0.5-1.0 ml/kg•h [6]) [7]. Yet, existing resuscitation protocols exhibit a large variability in treatment efficacy, due to many factors such as the inter-patient differences in response to burn injury and resuscitation as well as the incomplete knowledge of the pathophysiology underlying burn injury, which collectively complicates therapeutic decision-making. Hence, there is an ongoing effort to optimize burn resuscitation protocols [8]- [10]. However, due to the ongoing limitations, today's burn resuscitation often starts with an established burn resuscitation protocol but is subsequently titrated in an ad-hoc fashion to the physiological responses (including UO) in an individual patient to hopefully optimize the therapeutic efficacy. Unfortunately, current treatments tend to over-resuscitate the patients [11], which would exacerbate edema and expose the patients to an elevated risk of side effects, e.g., pulmonary edema, limb and abdominal compartment syndrome, necrosis, and death, due to the accumulation of resuscitation fluid (known as "fluid creep") [12], [13]. Hence, burn resuscitation regimens must be optimized based on complete knowledge of burn pathophysiology in order to best maintain organ functions in burn patients while minimizing adverse complications.
Development and optimization of burn resuscitation protocols is challenging, because (i) burn injury is less common than other widespread injuries (i.e., protocol optimization is hampered by the small number of burn injury cases) and (ii) it is unethical to test a new treatment protocol with unproven efficacy and safety in critically ill burn patients. From this standpoint, a mathematical model that can replicate physiological responses to burn injury and resuscitation in diverse burn patients can serve as an attractive basis to (i) perform non-clinical testing of emerging burn resuscitation protocols as medical digital twin [14], [15] and (ii) expand knowledge on burn pathophysiology [16]. However, existing mathematical models are not yet ideal for such purposes for at least three reasons. First, some mathematical models cannot be applied to test many existing burn resuscitation protocols based on UO, simply because they were not developed to predict UO responses to burn injury and resuscitation [17]- [21]. Second, most existing mathematical models have not been extensively proven in humans. In fact, the ability of existing mathematical models to replicate UO and other physiological responses to burn injury and resuscitation was validated in a prohibitively small number of patients [22]- [27] or only at the group level [17]- [21]. Third, some mathematical models (especially those reported early) [22]- [27] do not reflect up-to-date knowledge of burn-related physiology and pathophysiology gained in recent experimental and clinical investigations, regarding in particular the burn-induced perturbations in volume kinetics, kidney function, lymphatic flow, and tissue pressure-volume relationships [2], [28]- [30]. Hence, closing these gaps may lead us to an enhanced mathematical model of burn injury and resuscitation ideally suited to the development and testing of emerging burn resuscitation protocols and algorithms as well as to the expansion of our knowledge of burn pathophysiology.
In this paper, we intend to develop, extensively validate, and analyze a mathematical model capable of replicating volume kinetic and kidney function responses to burn injury and resuscitation in burn patients. By leveraging clinical datasets collected from 233 real burn patients receiving resuscitation, we developed a mathematical model suited to computer-aided in-human burn resuscitation trial and knowledge expansion, by expanding our prior work and utilizing systematic parametric sensitivity analysis and regularization. We investigated the validity of the mathematical model by testing its physiological plausibility in a dedicated test dataset. Using the validated mathematical model, we examined possible mechanisms responsible for the cohort-dependent differences in burn pathophysiology by comparing the mathematical models fitted exclusively to younger versus older patients, female versus male patients, and patients with versus without inhalational injury. To the best of our knowledge, our mathematical model may be the first mathematical model extensively validated for use as digital twin of real burn patients. This paper is organized as follows. Section II provides an overview of the mathematical model as well as the details of the clinical dataset and data analysis. Section III presents the results. Section IV is devoted to the discussion of major findings. Section V concludes the paper with future directions.
The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515

Mathematical Model
We continued to develop a mathematical model capable of predicting volume kinetic and kidney function responses to burn injury and resuscitation developed in our prior work [31]. The mathematical model consists of (i) volume kinetics to replicate water volume and protein concentration dynamics in the intravascular and the tissue compartments ("Volume Kinetics" in Fig. 1(a); see Appendix A), (ii) kidney functions to replicate UO response to changes in intravascular water and protein volumes ("Kidneys" in Fig.  1(a); see Appendix B), and (iii) transient perturbations in volume kinetics induced by burn as a chain of biochemical, molecular, and mechanical events ("Burn Perturbations to Volume Kinetics" in Fig. 1(b); see Appendix C).  Fig.  1(a) shows normal state, where the water (blue background) and protein (pink dots) balance between the intravascular compartment ("Plasma") and the tissue is preserved by the capillary filtration through the capillary pores (thin rectangles in the capillary wall) and lymph flow. The plasma is filtered by the renal plasma flow into the kidney by way of the renal regulatory mechanisms including the glomerulotubular balance and the Antidiuretic hormone. Upon the onset of burn injury, the water and protein volumes in the plasma decrease, while tissue compartment is divided into intact and burnt tissue compartments associated with water and protein content higher than normal state ("Burnt Tissue" and "Intact Tissue") ( Fig.1(b)). J: water flow. Q: albumin flow. Subscripts: C (capillary filtration); L (lymph flow); F (fluid infusion); U (UO); RPF (renal plasma flow); GFR (glomerular filtration rate); RR,GT (reabsorption rate by glomerulotubular balance); RR,ADH (reabsorption rate by Antidiuretic hormone); EX (exudation); EV (evaporation); PD (protein denaturation); BT (burnt tissues); IT (intact tissues). The yellow circles represent the perturbations that lead to the redistribution of water and protein. a: Partial destruction of capillaries in burnt tissues (shown as occluded capillary pores). b: Denaturation of protein in burnt tissues. c Transient negative hydrostatic pressure in burnt tissues, which draws water into the burnt tissues. d: Increased dermal fluid loss. e: Time-varying changes in capillary filtration and albumin permeability (shown as enlarged capillary pores). [f]: Vasodilation, which pushes water out of plasma. "X" ̅ represents "X" in normal (pre-burn) state.
The volume kinetics was represented using a three-compartmental model with intravascular space (including arterial and venous vessels, shown as "Plasma" in Fig. 1), intact tissues, and burnt tissues as separate compartments. These compartments describe the volumes of water ( , , ; whose dynamics is governed by Eq. (4)) and protein therein ( , , ; whose dynamics is governed by Eq. (5)) as functions of time. It was assumed that albumin serves as the surrogate of all the protein contributing to capillary filtration and colloid oncotic balance [32]. As shown in Fig. 1, the change in these volumes are dictated by (i) inter-compartmental flows of the lymphatic drainage (water and protein ; Eq. (8)-(9)) and capillary filtration ( and ; Eq. (6) and Eq. (10)); (ii) external inputs representing the gain via burn resuscitation ( and ); and (iii) outputs representing the kidneys' net filtration and reabsorption of renal plasma flow (i.e., UO ; Eq. (15) and Eq. (26)), as well as burn-induced evaporation ( ; Eq. (34)) and exudation ( and ; Eq. (35)). The kidney function was represented by a lumped-parameter model developed in our prior work [31], which includes a hybrid combination of first-principles and phenomenological elements that describe UO control governed by the kidneys, including the glomerular filtration rate (GFR; in Fig. 1 and Eq. (16)-(18)) modulated by the Starling forces due to the change in the intravascular water and protein volumes (i.e., plasma volume), the reabsorption by the glomerulotubular balance ( , in Fig. 1 and Eq.  25)). The transient perturbations in volume kinetics and kidney function triggered by burn injury were represented by an array of time-varying phenomenological models acting on various parameters and variables in the mathematical model of volume kinetics to replicate local and systemic pathophysiological changes caused by burn injury, including (i) partial destruction of capillaries in burnt tissues (a in Fig. 1; Eq. (28)-(29)), (ii) denaturation of protein in burnt tissues (b in Fig. 1; Eq. (31)), (iii) transient negative hydrostatic pressure in burnt tissues (c in Fig. 1; Eq. (32)), (iv) increased dermal fluid loss (d in Fig. 1; Eq. (34)-(35)), (v) time-dependent changes in capillary filtration and albumin permeability (e in Fig. 1; Eq. (28)-(30)), and (vi) vasodilation (f in Fig. 1; Eq. 36)). The disruption of the kidney function was the consequence of the perturbations occurring in volume kinetics as governed by Eq. (16)-(26). The equations associated with the mathematical model is summarized in Appendix.

Clinical Dataset
The clinical dataset used in this paper was furnished from two sources. The first source included 207 burn patients admitted to a burn intensive care unit (ICU) in December 2007-June 2009 [33]. These patients were treated with the aid of a clinical decision support system capable of recommending the hour-by-hour dose of lactated ringer (LR) to maintain UO at a target range of 30-50 ml/hr [10]. The care providers had the authority to override the recommendation. The dataset included hourly UO and LR dose as well as demographics including age, gender, and weight, total burn surface area (TBSA), the presence of inhalation injury, and the time of arrival. The second source included 53 burn patients. 29 patients were treated with the aid of the same clinical decision support system, while 24 patients were treated with the contemporary resuscitation protocols. The dataset included hourly UO and LR dose as well as demographics including age and weight, TBSA, and the time of arrival (gender and presence of inhalation injury were not known). Collectively, age, weight, and TBSA of the patients in the dataset were 47±18 years, 87±22 kg, and 40±18%, respectively. The overall mortality rate of the patients was 30%. In the first source, 77% of the patients were male and 11% of the patients were associated with inhalation injury.
The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 We randomly split the clinical dataset into training dataset to enhance and optimize the mathematical model (N=120) and test dataset (N=113) to validate the optimized mathematical model after excluding 27 burn injury patients associated with prohibitively small number of UO measurements (≤10). The demographic and injury severity of the burn patients in the training and test datasets were comparable (age: 45±19 years versus 49±18 years; weight: 85±18 kg versus 86±22 kg; TBSA 41.5±17.6% versus 38±18). The average hourly UO measurements in the training and test datasets were 23±2 samples and 20±4 samples, respectively.

Analysis
We conducted the analysis of the clinical dataset to continue to develop, optimize and validate the mathematical model to enable its in-human application (Fig. 2). First, we continued to develop the mathematical model by (i) including TBSA and weight dependency as well as human-compatible parameters to make it globally applicable to burn injury patients associated with diverse demographic characteristics and injury severity ("Model Enhancement" in Fig. 2), and then (ii) optimizing the mathematical model, by systematically reducing the parameter dimension using the training dataset, to make it compatible with sparse clinical measurements ("Model Optimization" in Fig. 2). Second, we validated the optimized mathematical model using the test dataset in terms of its predictive capability and physiological plausibility ("Model Analysis" in Fig. 2). Third, we scrutinized the mathematical models determined specifically for various categorical patient cohorts to gain insights on meaningful pathophysiological characteristics in these categorical patient cohorts ("Model Analysis" in Fig. 2). Details regarding the continued development, optimization, and analysis of the mathematical model are given in II.C.1, II.C.2, and II.C.3 below and are summarized in Table II.

Mathematical Model Enhancement
The clinical dataset presents major challenges in estimating all the parameters in the mathematical model on an individual patient basis. First, the burn patients in the clinical dataset exhibit large variability in the demographic characteristics as well as in the severity of burn injury (TBSA ranging between 16% and 100%), both of which increase the inter-individual variability in physiological responses to burn injury and resuscitation. Second, the number of parameters in the mathematical model is excessively large relative to the available measurements (i.e., hourly resuscitation dose and UO are the only measurements available to characterize these burn patients). To address these challenges and seamlessly apply the mathematical model to real burn patients, we leveraged the training dataset to further develop the mathematical model by (i) extending it to accommodate the variability in weight, TBSA, and species as well as (ii) systematically reducing the number of parameters to be estimated using the clinical dataset.
First, we categorized the parameters in the mathematical model into subject-invariant and subject-specific parameters. We categorized as subject-invariant parameters (i) those whose values appear consistent in multiple prior literatures (mostly associated with extensive properties and first principles components in the iii) The mathematical model was externally validated using testing dataset consisting of 113 burn patients. The mathematical model was analyzed to garner insights on the pathophysiological differences depending on age, gender, and inhalation injury. mathematical model, e.g., nominal water volume and albumin content in the intravascular and tissue compartments and the hydrostatic pressure in Bowman's capsule) and (ii) those whose values must be selected to yield mechanistically relevant physiological responses (e.g., parameters associated with the tissue compliance model, which must be chosen to result in physically relevant tissue hydrostatic pressure for a range of tissue volumes). The values of these subject-invariant parameters were mostly determined based on the existing literature (see Table AI for specific literatures we leveraged to determine these values). We categorized as subject-specific parameters (i) those whose values can exhibit large variability across burn injury patients (e.g., parameters pertaining to burn-induced perturbations and nominal glomerular filtration coefficient), (ii) those whose values have rarely been reported in the literature (e.g., capillary elastance and nominal lymphatic drainage rate), and (iii) those whose values are inherently unknown (e.g., parameters associated with phenomenological elements in the mathematical model). After all, a total of 58 parameters were categorized into 34 subject-invariant ("I" in Table AI) and 24 subject-specific ("S" and "SS" Table AI) parameters.
Second, we improved the mathematical model to accommodate the variability in weight and injury severity in the clinical dataset as well as to increase its suitability to real burn patients. To incorporate the weight dependence into pertinent parameters, we employed a linear allometric relationship by making them linear functions of weight, so that they assume typical values reported in the literature in case of a reference man (70 kg). These parameters include extensive parameters such as the water and protein volumes in the intravascular ( ̅ and ̅ in Table AI), intact tissue (V ̅ IT and ̅ in Table AI), and burnt tissue ( ̅ and ̅ in Table AI) compartments, capillary filtration rate ( ̅ in Table AI), and lymphatic drainage ( ̅ in Table AI) to list a few. One exception to the linear allometric relationship was the total body surface area ( in Table AI), which was made a function of weight through the Haycock formula (Eq. (34c)) and the weight-height relationship reported in the literature [34]. To incorporate the TBSA dependence into pertinent parameters, we (i) made the extensive parameters associated with the burnt tissue compartment functions of TBSA ( in Table AI) and (ii) expanded the plausible ranges of subject-specific parameter values associated with burn-induced pathophysiological responses so that the estimated parameter values avoid saturation at the pre-specified upper and lower bounds. These parameters include those representing the intensity of the inflammatory responses induced by burn injury such as the maximum increase in the capillary pore radius to albumin radius ratio (i.e. pore ratio) in the intact ( in Table AI) and burnt ( in Table AI) tissue capillary bed, the maximum drop in the burnt tissue hydrostatic pressure ( in Table AI), and the maximum increase in the capillary hydrostatic pressure ( in Table AI) to list a few. To make the mathematical model (which was initially developed based on the dataset collected from animals in our prior work [31]) more compatible with human burn patients, we refined the values of a number of parameters that are inherently different between animals and humans according to the literature (see Table AI for references). These parameters include nominal albumin concentration in the intravascular ([ ̅ ] in Table  AI) Table AI), and burnt tissue ([ ̅ ] in Table AI) compartments, colloid oncotic pressure constant ( in Table AI), nominal capillary hydrostatic pressure ( ̅ in Table AI), and the total body surface area ( in Table AI).

Mathematical Model Optimization, Training, and Validation
Using the training dataset, we optimized the mathematical model for in-human use by reducing the number of the subject-specific parameters that must be estimated on an individual patient basis. As described above, we down-selected 24 subject-specific parameters in the mathematical model that must be estimated using the hourly UO and LR dose measurements. Noting that the information content in the hourly UO and LR dose measurements may not be sufficient to robustly estimate all these 24 parameters, we capitalized on the training dataset to split the subject-specific parameters into parameters sensitive versus insensitive to the LR dose-UO input-output relationship. Then, we estimated the sensitive subject-specific parameters on an individual patient basis while fixing the insensitive parameters (together with the 34 subject-invariant The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 parameters) at their typical (i.e., group average) values. First, we determined the typical values of all the 24 subject-specific parameters by fitting the mathematical model to the LR dose-UO measurements pertaining to all the patients in the training dataset based on the pooled approach [35]. This task was accomplished by solving the following optimization problem using a multi-start gradient descent method ("globalsearch" in conjunction with "fmincon") in MATLAB (MathWorks, Natick, MA): where ̅ is the vector of typical values of subject-specific parameters, is the vector of 24 subject-specific parameters (i.e., a vector containing the 24 subject-specific parameters in Table AI), is the number of subjects, is the number of UO measurements for the subject during the treatment, ( ) is the value of UO associated with the subject measured at time , ( , ) is the value of UO associated with the subject at time predicted by the mathematical model for a given , and is the normalization factor for UO associated with the subject , which is defined as the range of multiplied by . Second, we classified the subject-specific parameters into sensitive and insensitive groups by quantifying and comparing the degree of inter-individual variability associated with all the subject-specific parameters. This task was accomplished by solving the following optimization problem for fitting with regularization [36] on an individual patient basis using a multi-start gradient descent method ("globalsearch" in conjunction with "fmincon") in MATLAB (MathWorks, Natick, MA): where is the vector of 24 subject-specific parameters associated with the subject , is the regularization weight, and is the normalization factor for the -th element ( ) of , which was defined so that all the elements in are homogeneously ranged approximately between 0 and 1. The regularization discourages the parameters from deviating from their respective typical values unless really needed to achieve superior goodness of fit. Hence, the subset of subject-specific parameters exhibiting deviations from the typical values in many subjects may be viewed as subject-specific parameters sensitive to the LR dose-UO input-output relationship. In this paper, we selected sensitive subject-specific parameters as those whose deviations exceeded a threshold value when averaged across all the 120 patients in the training dataset. Third, we ascertained the ability of the mathematical model (with the chosen sensitive subject-specific parameters) to faithfully replicate the UO responses to the LR dose in the training dataset, as well as its physiological plausibility. To this aim, we estimated the sensitive subject-specific parameters by solving the following optimization problem on an individual patient basis using a multi-start gradient descent method ("globalsearch" in conjunction with "fmincon") in MATLAB (MathWorks, Natick, MA) while fixing the remaining (insensitive subject-specific and subject-invariant) parameters at the respective typical values: where ̆ is the vector of sensitive subject-specific parameters selected by solving Eq. (2) (i.e., it is a subset of ), and ̆ is ̆ estimated for the subject . Then, we examined the faithfulness of the mathematical model in terms of (i) normalized mean absolute error (NMAE) [31], (ii) correlation coefficient, and (iii) UO range agreement, all between measured versus model-replicated UO on an individual patient basis, and (iv) Bland-Altman statistics between all measured versus model-replicated UO. We computed the UO The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 range-based agreement by specifying UO ranges of interest and then for each range computing the percentage of actual UO in the range whose model-predicted UO also resides in the same range. In addition, we examined the physiological plausibility of the mathematical model in terms of (i) both typical and subject-specific model parameter values (e.g., by comparing them with the values reported in the existing literature) as well as (ii) the plausibility of the volume kinetic and kidney function responses predicted by the mathematical model equipped with typical parameter values. We repeated the above procedure to optimize the mathematical model (i.e., Eq. (1)-(3)) so that it can yield minimal number of sensitive subject-specific parameters and adequate faithfulness and physiological plausibility.
Using the test dataset, we externally validated the faithfulness and physiological plausibility of the optimized mathematical model on an individual patient basis, in terms of the same metrics used above.

Mathematical Model Analysis
In addition to the optimization and validation of the mathematical model for in-human application (II.C.1-II.C.2), we also sought to garner in-depth insights and expand the knowledge on burn pathophysiology using the mathematical model. In particular, existing literature suggests that patients who are older [37], female [38]- [40], and associated with inhalation injury [37], [41], as well as those who are associated with severe burn injury [37] and/or receive delayed treatment [42] have a higher risk of mortality. The mathematical model already incorporates TBSA and arrival time post-burn, thereby allowing it to predict more severe responses to burn injury associated with large TBSA and delayed resuscitation treatments. However, it does not explicitly account for the effect of age, gender, and inhalation injury.
To investigate if the mathematical model can elucidate the age-, gender-, and inhalation injury-dependent differences in the burn physiological and pathophysiological mechanisms, we fitted the (optimized and validated) mathematical model separately to (i) younger versus older patients, (ii) female versus male patients, and (iii) patients with versus without inhalation injury. We used the patients in the training dataset (N=120), since they were associated with consistent treatment durations (i.e., 24 hours monitoring in most patients) compared to the test dataset. We excluded 16 patients since they did not have gender specification. We defined older patients as those with age above the median age of the 104 patients (45 years), and younger patients otherwise. Using these 104 patients in the training dataset, we built the group-average mathematical models associated with younger (N=52) versus older (N=52) patients, (ii) female (N=22) versus male (N=82) patients, and (iii) patients with (N=11) versus without (N=93) inhalation injury, all by solving a hybrid of the optimization problems in Eq. (1) and Eq. (3) (specifically, solving Eq. (1) only with respect to sensitive subject-specific parameters rather than all the sensitive parameters) based on the dataset associated with the specific patient groups. Then, we examined if the model parameter values for the two groups in each of the three categories (age, gender, and inhalation injury) exhibited meaningful contrasting differences that provide clinically important physiological insights.

Results
The iterative optimization of the mathematical model using the training dataset resulted in a mathematical model with seven sensitive subject-specific parameters in total, including the nominal capillary pore radius to albumin radius ratio ( ̅), the maximum increase in the pore ratio in the intact ( ) and the burnt ( ) tissue capillary bed, the maximum increase in the capillary hydrostatic pressure ( ), the slow decay rate associated with the increase in the capillary hydrostatic pressure ( 1, ), the tubule-glomerular feedback sensitivity ( ), and the nominal water reabsorption rate in the collecting ducts ( ̅ , ). Table I summarizes NMAE, correlation coefficient, UO range-based agreement pertaining to < 30, 30< <50, and >50 ml/h, or <0.5, 0.5< <1, and >1 ml/kg•hr, and the Bland-Altman statistics (i.e., the limits of agreement), all The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 associated with the optimized mathematical model. Fig. 3 presents examples of actual versus model-predicted UO responses associated with eight patients with various burn injury severity in the test dataset (details discussed in IV.A). Fig. 4 presents volume kinetic and kidney function responses to burn injury and burn resuscitation predicted by the group-average mathematical model in response to group-average burn resuscitation LR dose (details discussed in IV.A). Fig. 5 shows the weight-normalized PV, intravascular water gain (LR dose) and loss (capillary filtration in excess of lymphatic flow), and the burn resuscitation effectiveness (defined as the weight-normalized intravascular water gain rate (i.e., LR dose minus capillary filtration in excess of lymphatic flow) divided by the weight-normalized LR dose) throughout the 24-hour treatment period as predicted by the group-average mathematical model (details discussed in IV.B). Table II summarizes (i) demographics, (ii) statistical characteristics of fluid retention and UO relative to its treatment target range (30-50 ml/hr), and (iii) group-average model parameter values related to burn-induced inflammatory perturbations, all associated with the two patient groups in the three categories.

Discussion
Developing treatment strategies and expanding knowledge associated with burn injury present formidable challenges due to its complex pathophysiology, large inter-patient variability, and its less common incidence compared to other widespread injuries despite its devastating impact on the mortality and the quality of life. High-fidelity mathematical models capable of replicating volume kinetic and kidney function responses to burn injury and resuscitation has the potential to advance both treatment development and knowledge expansion aspects of burn resuscitation. Regardless, to the best of our knowledge, no mathematical model exists that has been developed and extensively validated using clinical datasets from real burn patients. In this paper, we present our continued development, extensive in-human validation, and analysis of a mathematical model for the study of burn injury and resuscitation, which is equipped with contemporary knowledge on the burn-related physiology and pathophysiology.

In-Human Credibility
The enhanced/optimized mathematical model based on the training dataset exhibited adequate predictive capability for UO response to burn injury and resuscitation in both training and test datasets (Table I and Fig.  3). In particular, the mathematical model worked equally well in both training and test datasets, both in terms of average statistics and robustness (e.g., NMAE was 15% with small IQR of 6%; Table I). Further, it could capture the physiological differences in burn patients across diverse TBSA range, including those associated with comparable weight ((a) versus (e) and (c) versus (g) in Fig. 3 Fig. 3). In addition, the mathematical model showed a high degree of UO range-based agreement (>90% (when weight-normalized) and >78% (when not weight-normalized) of model-predicted UO resided in the same range to actual UO; Table I). Noting that existing burn resuscitation protocols determine the hourly resuscitation dose based on the range of UO, the results suggest that the mathematical model may serve as a valuable platform for non-clinical testing of burn resuscitation protocols and algorithms. In addition to UO, the mathematical model was able to predict the overall volume kinetic and kidney function responses to burn injury and resuscitation in a realistic way: the behaviors of the internal volume kinetic and kidney variables were consistent with the contemporary knowledge on burn pathophysiology as well as findings from recent studies (Fig. 4).

) and those associated with distinct weight ((b) versus (f) and (d) versus (h) in
The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515   Specifically, the group-average mathematical model predicted that (i) plasma volume and UO showed an anticipated trend of initial decline upon the onset of burn injury and subsequent recovery with the initiation of burn resuscitation and later with the return of resuscitation fluid leaked into the tissues back to the blood ( Fig.  4(b), Fig. 4(c), and Fig. 4(e)) [43], [44]; (ii) burnt and intact tissue volumes increased up to nearly twice their initial values, peaking and starting to decay approximately at 24 hours post-burn ( Fig. 4(b)) [45]- [47]; (iii) plasma albumin was transported into burnt and intact tissues due to the perturbations in albumin reflection and permeability-surface area coefficients (Fig. 4(d) and Fig. 4(k)) triggered by burn-induced increase in the capillary pore size that decreased the capillary pore radius ratio in both burnt and intact tissue (Fig. 4(i)) [48]; (iv) GFR increased just a few hours post-burn even before plasma volume (Fig. 4(f)) [49]; and (v) sodium concentration decreased after burn injury and resuscitation (Fig. 4(h)) [50].
Importantly, the mathematical model could predict UO as well as physiologically plausible volume kinetic and kidney function responses once physiologically acceptable values were assigned to its parameters. In fact, the majority of the parameters equipped with physiological implications assumed values comparable to typical values and/or those reported in the literature both on the individual and population-average basis (Table AI).
In sum, we demonstrated that the mathematical model can faithfully replicate the volume kinetic and kidney function responses in a wide range of burn patients, both in terms of the adequacy of the model-predicted responses and the plausibility of the model parameter values.

Insights on Burn Resuscitation Effectiveness
One strength of the mathematical model presented in this paper is its ability to replicate overall responses of a burn patient to injury and resuscitation, including those that cannot be clinically measured. Exploiting this advantage, we sought to garner insights on the effectiveness of burn resuscitation in a typical patient subject to burn injury. It is known that the homeostasis in volume kinetics is severely disrupted after burn injury due to the activation of multiple inflammatory mediators, which in turn causes a large portion of the resuscitation fluid to leak out of the intravascular compartment via capillary filtration. Although this leakage is partially recovered by the increase in lymphatic flow, >50% of the resuscitation fluid can leak out of the intravascular compartment in the initial hours post-burn in extensive burn injury [51], [52]. In this regard, burn resuscitation effectiveness represents the portion of the resuscitation fluid actually used to expand plasma volume. Based on the investigation and interpretation of >50 physiological variables including those presented in Fig. 4 and Fig. 5, we could garner the following insights on the important physiological and pathophysiological mechanisms responsible for the effectiveness of burn resuscitation during the initial 24 hours post-burn. Initially, there is a large and fast fluid shift from the blood to the intact and burnt tissues immediately after burn injury for up to one hour, leading to a large decrease in the plasma volume (Fig. 4(j) The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 and Fig. 5(a)). Our mathematical model suggests that major mechanisms responsible for this initial loss of plasma volume may be negative hydrostatic pressure, protein denaturation in burnt tissues, and systemic increase in the capillary hydrostatic pressure.
After this initial phase, a decrease in the capillary filtration and the self-regulation of plasma volume occurs for up to one hour. Our mathematical model suggests that major mechanisms responsible for this phase may be the reduction in the plasma volume and the resulting decrease in the capillary hydrostatic pressure, the recovery of hydrostatic pressure in the burnt tissues, and the increase in the lymphatic flow (Fig. 4(l)). The effectiveness of burn resuscitation during this phase is very high (>100%), meaning that plasma volume is expanded based on almost all the resuscitation fluid as well as the fluid returning from the edematous (burnt and intact) tissues (Fig. 5(c)). Subsequently, burn resuscitation effectiveness is deteriorated quickly as the capillary filtration of water and protein increases again due to the opened capillary pores (Fig. 4(i)) and the increase in the plasma volume. Our mathematical model suggests that major mechanisms responsible for this phase may include the increase in the protein concentration in both burnt and intact tissues as well as the hypoproteinemia in the blood (Fig. 4(d)), which altogether increase the osmotic pressure gradient toward tissues and promote capillary filtration of both water and protein (thereby forming a vicious circle). Burn resuscitation effectiveness reaches its minimum level of 2%-15% at 10-15 hours after the initiation of treatment, which is in close agreement with the literature suggesting maximal edema formation in this period post-burn [45]. Finally, burn resuscitation effectiveness increases back to approximately 40% at 24 hours after the initiation of treatment. Our mathematical model suggests that major mechanisms responsible for this recovery may include the recovery of lymphatic flow to return excessive water and protein to the blood as well as the gradual decrease in the capillary pore size, which altogether decreases the fluid extravasation rate.

Cohort-Dependent Differences in Burn Pathophysiology
The analysis of datasets associated with various categorical patient cohorts (with respect to age, gender, and inhalation injury) provided meaningful insights on the cohort-dependent differences in burn physiology and pathophysiology (Table II). To begin with, the mathematical model was able to replicate UO response to burn injury and resuscitation associated with all the categorical patient cohorts (younger versus older, female versus male, and patients with and without inhalation injury).
First, between younger versus older patients, the latter had much higher mortality rate and higher portion of UO responses below the target therapeutic range than the former despite its smaller group-average TBSA and higher level of weight-normalized LR dose, which leads to higher fluid retention in the latter (Table II, first two columns). Comparing the mathematical models fitted to younger versus older burn patients, the latter was associated with higher inflammation factors (including the larger increase in capillary pore size in the burnt tissues (M α BT ) and the capillary hydrostatic pressure (M P C )). Second, between patients with versus without inhalation injury, the former likewise had >1.3 times higher mortality rate and higher portion of UO responses below the target therapeutic range than the latter (Table II, last two columns). Comparing the mathematical models fitted to burn patients with versus without inhalation injury, the former exhibited higher inflammation factors (including the larger increase in capillary pore size ( and ) and capillary hydrostatic pressure ( )) similarly to the older patient cohort. Existing literature shows the possible association between fluid retention and edema versus mortality and complication rates in burn patients [12]. In addition, both literature and our dataset indicate higher fluid retention and edema as well as higher mortality rates in older patients and patients with inhalation injury [53], [54]. From this standpoint, our mathematical model analysis predict that higher inflammation may be a key contributing factor in increasing the mortality risk of older patients and patients with inhalation injury via an increase in the fluid retention.
The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 Our prediction is in fact consistent with the contemporary knowledge in the literature identifying inflammation as an important mediator of increased fluid retention and edema with the increased mortality rate in older patients [55], [56] and patients with inhalation injury [57], [58], although other causes can play a role (e.g., degraded cardiovascular efficiency in elderly burn patients [59]). Although our mathematical model analysis reveals possible mechanisms responsible for higher mortality rate, the exact cause is yet to be clearly elucidated. Regardless, lower burn resuscitation effectiveness in older patients and patients with inhalation injury relative to younger patients and patients without inhalation injury remains true, and our mathematical model was able to replicate the age-and inhalation injury-dependent differences in burn resuscitation effectiveness. Hence, our mathematical model may serve as an effective basis to develop and validate burn resuscitation protocols and algorithms suited to these categorical patient cohorts.
Third, between female and male patients, neither the dataset nor the mathematical model showed any meaningful difference in terms of inflammation and fluid retention (Table II, two middle columns). This contrasts against some literature identifying the female gender as a mediator of mortality risk associated with burn injury, as confirmed by our dataset (1.3 times the mortality rate in males). Hence, our mathematical model analysis suggests that higher mortality risk in female burn patients may be attributed to factors other than an increase in the inflammation and the corresponding increase in the fluid retention, especially those not manifested in the initial 24 hours post-burn. In fact, a prior study performed on mice showed that the difference in the inflammatory responses in female and male subjects was not clear until 6 days post-burn [60]. The exact mechanisms responsible for the gender difference in burn-induced mortality risk are still unknown and must be unveiled. Table II: Demographics, characteristics of fluid resuscitation, fluid retention, and urinary output (UO) relative to its treatment target range (30-50 ml/hr), and group-average model parameter values related to burn-induced inflammatory perturbations, all associated with the two patient groups in age (younger versus older patients), gender (female versus male patients), and inhalation injury (patients with versus without inhalation injury) categories. Fluid retention is computed as the total resuscitation fluid (LR) volume minus the total UO during the 24 hours of treatment. Fluid resuscitation and fluid retention are shown in mean±SD, while UO is shown in median (IQR).
, , and are group-average values of burn-induced inflammatory perturbation parameters in the mathematical model. The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515

Conclusions
The main outcome of this paper is a physiologically plausible mathematical model capable of replicating volume kinetic and kidney function responses to burn injury and resuscitation suited to in-human application. To the best of our knowledge, the mathematical model presented in this paper may be the first of its kind developed and extensively validated using large clinical datasets from real burn patients. We anticipate that the mathematical model may provide an attractive platform to conduct non-clinical testing of burn resuscitation protocols and test new hypotheses on burn pathophysiology. Future effort must be exerted to investigate the potential of the mathematical model as medical digital twin for disciplined development and rigorous stress testing of emerging burn resuscitation algorithms and as a cornerstone to expand our understanding of burns.

Appendix
This section summarizes the mathematical model equations and parameter values. Full complete details on the mathematical model can be found in our prior work [31] and the references therein.

Volume Kinetics
Volume kinetics is represented by a classical multi-compartmental model consisting of vasculature, intact tissues, and burnt tissues as separate compartments. The water (Eq. (4)) and albumin (Eq. (5)) balances in each compartment are formulated based on the mass conservation principle: where is water volume, is albumin content, is water flow, and is albumin flow. The subscripts , , , , , , and denote capillary filtration, lymphatic flow, fluid infusion, UO, evaporation, exudation, and protein influx due to burn-induced denaturation (see Appendix C), respectively, while the subscripts P, BT, IT denote plasma, burnt tissues, and intact tissues. The capillary filtration is expressed using the Starling equation: where ∈ { , }, and are capillary and tissue hydrostatic pressures, and are plasma and tissue colloid oncotic pressures, and , and are capillary filtration and albumin reflection coefficients, respectively. [ ] and [ ] are plasma and tissue albumin concentrations, respectively.
is the colloid The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 oncotic pressure constant. The lymphatic water flow is expressed by a phenomenological model in the form of a sigmoidal curve [31]: where ̅ , is nominal lymphatic flow at nominal tissue hydrostatic pressure ̅ , and and are constants representing the inverse of the maximal degree of increase in the lymphatic flow and its sensitivity to the change in the tissue hydrostatic pressure. The lymphatic albumin flow depends both on the lymphatic water flow and the albumin concentration of the corresponding tissue: The albumin transport across the capillary wall is expressed based on the coupled diffusion-convection equation [61], where PS is permeability-surface area coefficient: A linear phenomenological relationship between P C and V P is assumed: where ̅ is nominal capillary hydrostatic pressure associated with the nominal PV ( ̅ ), and is capillary elastance. For the interstitial compartments, we use the nonlinear microscopic pressure-volume model developed by Øien and Wiig [62] and employ it in our macroscopic context: where , , ̂ , , and n are constant coefficients (Table AI), is a measure of hydration in the tissue at the microscopic level, and ̂ and ̌ its maximum and minimum values, and R is the radius of glycosaminoglycans, varying in response to . Eq. (11) relates the microscopic hydration to macroscopic hydration, where ̅ is nominal hydration level, defined as the ratio of the water volume in and its dry weight [30].

Kidney Function
The renal function is based on by a novel lumped-parameter, hybrid mechanistic-phenomenological model we developed in our previous work [31]. As blood circulates in the body, the renal plasma flow ( ) flows into the kidneys to be filtered. The glomerular filtration rate ( ) is the fraction of that is filtered into the Bowman's capsule to flow along the nephron tubules, the majority of which is reabsorbed back into the circulation, cleared from waste material. UO is given by the difference between glomerular filtration rate ( ) and reabsorption rate (RR; ): The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 is derived from the Starling forces between glomerular capillaries and Bowman's capsule as follows [31]: where is a constant coefficient relating the glomerular hydrostatic pressure ( ) to the glomerular resistance ( ) and , and is the glomerular filtration coefficient. The hydrostatic pressure in Bowman's capsule ( ) is assumed constant during burn injury and resuscitation [63]. Assuming that the albumin flowing into the kidney does not pass into the Bowman's capsule, we can write: (16.c) where [ ] is the colloid content of the glomeruli, and = is the filtration fraction. Since is naturally always less than 1, assuming the colloid oncotic pressure constant in the glomeruli is identical to other capillaries: is perturbed by the fluctuations in cardiac output, which in our mathematical model is assumed to be proportional to PV: where ̅ is nominal corresponding to nominal plasma volume ̅ . The perturbation in is strictly compensated for by the modulation of renal capillary resistance via tubulo-glomerular feedback (TGF). We express TGF as simple phenomenological dynamic modulation of which is a hypothetical variable representing glomerular resistance [31]: where = − ̅ is the deviation of from its nominal value ̅ enforced by TGF, and are time constant and sensitivity associated with TGF, and ̅ is the nominal value of . Combining Eq. (16)- (18), is expressed in terms of PV and plasma albumin content.
is assumed to be modulated by the glomerulotubular balance (GT, , ), and pure water reabsorption in the collecting ducts governed by the Antidiuretic hormone levels (ADH, The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 GT modulates proximal tubules so that approximately 80% of GFR is reabsorbed [64], [65]: ADH content is modulated by baroreceptor (inversely related to the changes in PV) and osmoreceptor (directly related to the changes in sodium concentration ([Na + ]) in the blood) signals [66], [67]. Typically, 27% of ADH is excreted without reabsorption via GFR, while 73% is metabolized in the liver proportionally to blood flow (which is assumed to vary with PV). The dynamics of the ADH content may then be expressed as follows [31]: where ADH is ADH content in the extracellular fluid, is nominal ADH secretion rate, and [ + ] are positive constant coefficients, and [ ] is ADH concentration in the extracellular fluid: where is maximum , and [ ] 50 is ADH concentration corresponding to = 1 2 . Based on the idealistic assumption that sodium content is conserved in the body (including the collecting ducts), and is diluted as the pure water reabsorption increases in the collecting ducts [67]: where [ + ] ̅̅̅̅̅̅̅̅ is nominal sodium concentration corresponding to ̅ , . Combining Eq. (15)- (25), is expressed in terms of plasma volume, plasma albumin content, and ADH concentration:

BURN-INDUCED PATHOPHYSIOLOGY IN VOLUME KINETICS AND KIDNEY FUNCTION
To express various perturbations induced by burn injury, a universal perturbation function below is used: where is maximum perturbation, and 1, and 2, are slow and fast time constants, all corresponding to a perturbation ∈ { , , , }, where and are pore radius ratios associated with burnt and The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 DOI: 10.1109/TBME.2021.3094515 intact tissues (see below for details). Eq. (27) describes an intensifying perturbation post-burn, which reaches its maximum at a varying time post-injury and then starts to decay. It can be incorporated into volume kinetics and kidney function models above to replicate burn-induced perturbations. Details follow.
The destruction of capillaries in burnt tissues is expressed as perturbations in capillary filtration coefficient , and permeability surface area coefficient . Based on the pore theory of trans-capillary exchange [68], [69], capillary filtration, albumin reflection, and permeability-surface area coefficients associated with burnt and intact tissues are expressed as functions of the pore ratios , ∈ { , }, defined as the ratio between albumin radius and capillary pore radius: where ̅ , , ̅̅̅̅ , and ̅ are nominal values of , , , and adjusted for the water fraction in burnt and intact tissues, weight, and capillary recruitment [21]: where ̅ and ̅̅̅̅ are nominal capillary filtration and permeability surface area coefficients in the absence of burn injury, = (2 ̅ − 1) (which shows that the parameters in Eq. (29) are extensive and depend on PV), and is the fluid volume ratio between skin and total interstitial compartment [21]. Burn-induced changes in these coefficients are expressed by formalizing the burn-induced changes in the capillary pore radius ratios using Eq. (27): The denaturation of protein in burnt tissues ( in Eq. (5b)) is expressed as a protein influx into burnt tissues [28], [48], [70]: where ̂ is protein influx immediately post-burn, decaying with a time constant .
The transient negative hydrostatic pressure in burnt tissues is expressed using Eq. (27): The overall hydrostatic pressure in burnt tissues is computed by combining Eq. (12) and Eq.
Evaporation and exudation (Eq. 4(b)-4(c)) are modeled based on the literature [71], [72]: The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 where 1, , 2, , , 1, , 2, , and are constant coefficients (Table AI), is the fraction of body surface area subject to burn, is total body surface area based on the Haycock formula, and is the ratio between albumin concentration in the exudate and albumin concentration in the burnt tissues [21].
The perturbation in capillary hydrostatic pressure is expressed using Eq.
The overall capillary hydrostatic pressure is computed by combining Eq. (11)   The final version of this paper has been accepted for publication in the IEEE journal of Transactions on Biomedical Engineering and is available for early access at https://ieeexplore.ieee.org/document/9478222 -