Mathematical modeling of the combined effects of thermal burn and local irradiation

TL;DR

This paper introduces a new mathematical model to study the combined effects of radiation and thermal burn injuries, improving predictions of their impact.

Contribution

A novel mechanistic model capturing synergistic effects of radiation and burn injury with limited data.

Findings

The model successfully reproduces inflammatory and fibroblast behaviors seen in thermal and combined injuries.

Sensitivity analysis and parameter sampling enabled effective use of limited data for modeling.

The framework allows for better exploration of inflammatory mechanisms in radiation combined injury.

Abstract

Radiation combined injury (RCI), resulting from ionizing radiation exposure accompanied by other injuries such as burn, laceration, or fracture, is associated with higher rates of mortality and severe effects. Current available models lack the ability to capture synergistic effects associated with RCI and/or rely on data-driven approaches that are limited in their predictive capabilities. To address this, we developed a mechanistic mathematical model for local radiation exposure combined with burn injury that captures the inflammatory response and early fibroblast activity associated with injury resolution. We utilized sensitivity analysis and parameter sampling methods to leverage limited data. The model was able to reproduce inflammatory and fibroblast behavior consistent with observed thermal injury and combined injury profiles. The formulation of a mechanistic model for combined injury adds increased modeling flexibility allowing for further exploration of the underlying inflammatory mechanisms and provides a framework for leveraging minimal data for improved predictive models of RCI.

Figures

Schematic diagram of locally irradiated superficial thermal burn, the state variables, and their relationships.The model consists of the thermal burn compartment where cell populations experience damage from thermal fluence and the surrounding tissue compartment where resident immune cell populations experience damage from a prompt radiation dose. The variables model the inflammatory and early proliferation phases of wound healing. State variables are represented using circles and squares. The solid black lines indicate transitions between different cellular states. Solid lines that feedback into the state variable indicate proliferation. Initiated cellular interactions are captured by dotted lines and remain locally concentrated. Blue and white colors are used to distinguish between upregulation (i.e., increasing in response to the cells themselves) and destruction/inhibition, respectively. The toxic triangles indicate cell transitions and interactions affected by ionizing radiation that have been incorporated into the model.

M1 macrophage ratio to total macrophages and damage variable with varying prompt radiation doses.The M1 ratio (panel A) denotes the total expression of M1 macrophages as a ratio of the total activated macrophage response such that 1 corresponds with 100% M1 activated macrophages and 0 corresponds with 0% M1 activated macrophages (100% M2 activated macrophages). Panel B shows the corresponding damage variable indicating damage resolution over time. Prolonged pro-inflammatory behavior has been widely evidenced to delay the wound healing process [48].

Plots of damage (Damtb) and debris (Debtb) for a superficial thermal burn with varying levels of radiation. Panel A shows Damtb and panel B shows Debtb. Superficial thermal burns alone should heal fully within 4 to 5 days. The introduction of ionizing radiation delays the immune response and results in prolonged healing times.

Plots of undamaged cells in the surrounding tissue supplied by the bloodstream.Panel A represents neutrophils, panel B represents macrophages, panel C represents lymphocytes, and panel D represents fibroblasts. Populations in the surrounding tissue are affected by both the initial cellular injury from ionizing radiation as well as reduced influx due to microvasculature damage in the surrounding tissue.

Plots of immune cells in the injury compartment.Panel A represents neutrophils, panel B represents M1 macrophages, panel C represents L1 lymphocytes, panel D represents M2 macrophages, and panel E represents L2 lymphocytes. Immune cell response in the injury is delayed due to ionizing radiation exposure.

Plots of fibroblasts (Ftb) and damage (Damtb) in the thermal injury.Ionizing radiation exposure alters fibroblast function resulting in reduced proliferative abilities and irregular collagen deposition, delaying damage resolution.

Plots of pathogen (Ptb) and damage (Dam_tb) at the thermal injury site.Since ionizing radiation is associated with risk of bacterial infiltration leading to infection, the model includes the capability to capture prolonged bacterial presence and delayed healing which could result in a non-healing wound and chronic inflammation.

Equations

Full text

Introduction

In the context of nuclear accidents, such as Chernobyl, radiation combined injury (RCI) describes a condition wherein radiation exposure is combined with additional injuries such as burns, wounds, and other trauma [1–3] Combined injuries are associated with greater tissue damage and an increased risk of adverse outcomes, such as lethality [4,5], secondary infection [6], and sepsis [7–9]. Moreover, synergistic effects have also been observed [4,9] where the combined effects are more severe than the sum of the individual injuries.

The pathophysiology of individual burn injuries has been adequately captured using mathematical models [10–12], which determine the resulting tissue damage and underlying cellular activity during injury resolution. Similarly, mathematical models of injury from ionizing radiation have also been developed [13–18]. These models typically address radiation at the cellular level and provide estimations for direct cell damage, myelopoietic recovery, and survival. However, a combined injury model that captures the synergistic pathophysiologic effects between a thermal burn combined with radiation is not currently available. This is especially crucial since local ionizing radiation exposure has been observed to delay the burn injury wound healing process [19–22] resulting in depleted circulating immune cells and additional damage to resident cell populations. Local ionizing radiation exposure has also been evidenced to damage dermal fibroblasts [23] and microvasculature [24], which provide additional causes of delay and dysregulation of the inflammatory response [25,26]. These effects of both the burn injury and local radiation exposure contribute synergistically to delayed resolution of injury and inflammatory imbalances.

Given the lack of mechanistic mathematical models, empirical models have been primarily used to estimate injury prevalence [27] and probability of survival [28]. However, empirical models are bound by their data and rarely maintain accuracy when extrapolated beyond recorded time points. Alternatively, mechanistic models can capture system dynamics that allow for exploration past the constraints of available data as well as the capability to estimate related effects such as medical outcomes associated with the injury.

For our model, we developed a system of ordinary differential equations (ODEs) which capture the immune response to radiation combined injury with burn. A variety of methods have been used to mechanistically model the immune response to injury including systems of ODEs [29–32], partial differential equations (PDEs) [33,34], and agent-based models [35–37]. Each method has advantages and disadvantages, mainly regarding biological complexity and computation time. PDEs and agent-based models, for example, are often chosen due to their ability to incorporate spatial gradients important in immune cell infiltration at the wound site. However, this added complexity can often result in increased computation time. For this reason, we have selected a system of ODEs that can sufficiently capture cellular population dynamics within the injury site and surrounding tissues while keeping computation time reasonable. This supports future integration into predictive software for modeling combined injury in nuclear incidents.

To model the effects of thermal burn combined with local radiation, we developed a time-dependent physiological model (TDPM), which includes immune cell infiltration, inflammation, and early fibroblast activity. We employed a combination of data, medical trends, and sensitivity analysis to generate plausible parameter sets and conducted numerical simulations to demonstrate the ability of the model to reproduce observed combined injury dynamics.

Methods

The mathematical model is a system of ODEs representing the various immune cell populations, damaged tissue and cellular debris, and foreign pathogen. A system of ODEs was chosen to capture nonlinearities and feedback loops present in inflammatory signaling and immune cell response while keeping computation time reasonable. ODEs, however, cannot account for spatial complexity and assume the elements in the environment are well-mixed and evenly distributed throughout the space. One way to resolve this is to use a compartmental model in which cells and cellular mediators can travel between compartments. This method is employed in the present model.

In the mathematical model for a locally irradiated thermal burn injury, we consider two compartments: 1) the thermal burn compartment with activated immune cells, fibroblasts, and foreign pathogen, and 2) a surrounding irradiated tissue compartment, which serves as an intermediate compartment between circulating immune cells in the bloodstream and their activation at the injury site. Cell populations and their associated cytokine and chemokine biomolecular signals are represented with compartmental ODEs. A model schematic of the variables and interactions can be viewed in Fig 1. Detailed descriptions of model variable descriptions and symbols are available in Table 1. Parameter descriptions along with numerical values used in model simulations are provided in S1 Table.

Table 1: State variables with descriptions.

Schematic diagram of locally irradiated superficial thermal burn, the state variables, and their relationships.The model consists of the thermal burn compartment where cell populations experience damage from thermal fluence and the surrounding tissue compartment where resident immune cell populations experience damage from a prompt radiation dose. The variables model the inflammatory and early proliferation phases of wound healing. State variables are represented using circles and squares. The solid black lines indicate transitions between different cellular states. Solid lines that feedback into the state variable indicate proliferation. Initiated cellular interactions are captured by dotted lines and remain locally concentrated. Blue and white colors are used to distinguish between upregulation (i.e., increasing in response to the cells themselves) and destruction/inhibition, respectively. The toxic triangles indicate cell transitions and interactions affected by ionizing radiation that have been incorporated into the model.

Model equations

The system of ODEs were developed using terms and model structure similar to those found in previous models of inflammation [29–31] and proliferation [32] with modifications made to incorporate radiation effects, which are discussed herein. Additional details regarding the development of specific model terms can be found in S2 File. A few auxiliary functions (Equations 1–3) are used throughout the system of equations which represent an inhibitory function, $[eqn]$ , a hill function, $[eqn]$ , and a dose-dependent saturating function between 0 and 1, $[eqn]$ , respectively. The inhibitory function and hill function are commonly used across models of inflammation [29–32].

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Damaged tissue and debris

Damaged tissue and cellular debris result from direct damage to cells from thermal fluence as well as collateral damage due to proteases and other tissue degrading molecules produced by neutrophils and M1 macrophages during the inflammatory phase [38–40]. The resulting cellular debris is then removed through phagocytosis by activated inflammatory cells (neutrophils and macrophages). Damaged tissue is then resolved through fibroblast proliferation which serves as the first step towards replacement of damaged tissue.

$$[eqn]$$ $$[eqn]$$

Although collagen is not explicitly modeled, damage resolution by activated fibroblasts captures early repair processes indicated by collagen deposition at the wound site. Collagen has been observed to be a good marker of tissue damage resolution in other mathematical models of wound healing [32]. Radiation has been evidenced to inhibit the deposition of collagen by fibroblasts [22,41,42], which is included as an inhibitory function in the repair term for the damage equation. The population of damaged fibroblasts in the surrounding tissue ( $[eqn]$ ) inhibits the removal of damage.

Neutrophils

Neutrophils are typically the first inflammatory cells to arrive at the wound site and are activated and recruited by damage-associated molecular patterns and other immune cells [40]. Neutrophils then work to phagocytize any debris and foreign matter.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Localized ionizing radiation exposure damages the microvasculature [22,26] resulting in delayed infiltration of recruited neutrophils and other immune cells [25]. Revascularization is not explicitly modeled, so this delay is represented by a dose-dependent inhibitory function that decays over time and reduces the baseline influx of neutrophils.

Monocytes and macrophages

Monocytes and macrophages are also recruited to the injury site to remove any debris and regulate the inflammatory response [40]. In our model, we use the term monocyte to represent the inactive form and macrophage to represent the activated form. Monocytes in the surrounding tissue are present at rest and have two states: undamaged or damaged, denoted by a $[eqn]$ or $[eqn]$ superscript, respectively. Macrophages in the injury site can be activated on a spectrum presenting a range of pro- and anti-inflammatory behavior [43], which are generally labeled as M1 classically-activated and M2 alternatively-activated respectively. Both are used in the model to represent these general phenotypes.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Undamaged monocytes in the surrounding tissue are recruited from the blood stream at a baseline constant rate, represented by the parameter $[eqn]$ in the first term of Equation 9. Similar to the equation defined for neutrophils in the surrounding tissue (Equation 6), this baseline rate is inhibited when radiation injury occurs in a dose-dependent manner.

For prompt gamma radiation doses above 2 Gray (Gy), macrophages are polarized towards the classically activated M1 phenotype [44]. Given the prolific evidence of increased pro-inflammatory cytokines following irradiation [45–47], it is reasonable to assume that this overexpression of pro-inflammatory cytokines could be the mechanism by which activation towards the M1 phenotype occurs. For this reason, the alternative activation is dependent upon damaged macrophage cells in the surrounding tissue which includes all associated cytokines, chemokines, and other signals. This allows for both a time- and dose-dependent inhibition of the activation process toward the M2 phenotype. This can be observed in Fig 2. Panel A shows the scaled M1 ratio to total macrophages in the thermal burn compartment while panel B displays the corresponding damage variable over time. The plots exhibit the increased and prolonged activation toward the M1 phenotype with increasing radiation dose as well as the delayed wound healing exhibited by prolonged inflammation and delayed resolution of damage. This mechanism is represented in the first term of Equation 12. Consistent with current research [44], the activation towards the M2 phenotype is only inhibited for $[eqn]$ Gy, otherwise, activation occurs at the rate $[eqn]$ .

M1 macrophage ratio to total macrophages and damage variable with varying prompt radiation doses.The M1 ratio (panel A) denotes the total expression of M1 macrophages as a ratio of the total activated macrophage response such that 1 corresponds with 100% M1 activated macrophages and 0 corresponds with 0% M1 activated macrophages (100% M2 activated macrophages). Panel B shows the corresponding damage variable indicating damage resolution over time. Prolonged pro-inflammatory behavior has been widely evidenced to delay the wound healing process [48].

T lymphocytes

T lymphocytes also infiltrate the injury during the inflammatory phase to aid in debris removal and inflammatory regulation, but typically at later time points than neutrophils and macrophages. Like monocytes, lymphocytes have inactive populations in the surrounding tissue and are highly radio sensitive [49,50]. Populations in the surrounding tissue will also be split between undamaged and damaged, denoted by a $[eqn]$ or $[eqn]$ superscript, respectively. Only T-cell lymphocytes are considered in the present model which include the pro-inflammatory phenotypes (e.g., $[eqn]$ T cells, TH1 cells, and T17 cells) and anti-inflammatory phenotypes (e.g., TH2 cells and regulatory T cells).

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Lymphocyte function in the present model consists of general regulation of the macrophage response through inflammatory mediators and other signaling biomolecules. At present, damaged lymphocytes in the surrounding tissue do not affect lymphocyte function or polarization, as was seen in the macrophage equations, but the population of damaged lymphocytes is included due to their importance in models of radiation.

Fibroblasts

The presence of fibroblasts in the injury site indicates the start of the proliferative phase which is marked by the deposition of new extracellular matrix to replace damaged tissue [40]. Fibroblasts in the surrounding tissue are recruited from the local connective tissues as well as through proliferation of resident populations.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Recruitment of fibroblasts from the surrounding tissue is inhibited like the inflammatory cells in a dose-dependent manner to account for the microvascular damage associated with local tissue injury from radiation. Fibroblast proliferation is also inhibited by radiation injury as this has been observed to reduce the proliferative abilities [23,51]. Similar to previous equations, we use the population of damaged fibroblasts to drive this inhibition.

Equation 22 has an analytical solution solved using separation of variables and simple integration. Since this equation only has constant parameter rates, the final solution is only dependent on time. The solution is represented in Equation 24 and was used throughout simulations in place of Equation 22.

$$[eqn]$$

Pathogen

Ionizing radiation is associated with an increased risk of infection regardless of a localized [52] or whole-body exposure [53–56]. This risk is observed to be further elevated for irradiated burns [7] so pathogen dynamics must be accounted for in the model. The pathogen equation was generated using a common equation structure for bacterial populations [30]. The equation includes logistic growth and removal by phagocytosis from the model’s inflammatory cells and other background responses (e.g., mast cells and natural killer cells).

$$[eqn]$$

Initial conditions

The model variables are initialized using a combination of thermal fluence, ionizing radiation dose, and steady state calculations. The cellular damage in the thermal burn injury is assumed to be a result of thermal fluence only and the damage in the surrounding tissue is assumed to be a result of radiation exposure only.

Surrounding tissue compartment.

The surrounding tissue compartment contains tissue resident immune cell populations and fibroblasts unaffected by the burn injury. Neutrophils have not been evidenced to have tissue resident populations [57], so this variable is initialized at zero. Monocytes, lymphocytes, and fibroblasts do have tissue resident populations at a homeostatic steady state within the system. The steady value for these variables is the ratio of the constant influx rate from the blood vessels to their decay rate. Thus, the steady state in the surrounding tissue for monocytes is represented by $[eqn]$ and the steady state for lymphocytes is represented by $[eqn]$ The fibroblast population will approach the ratio of the influx rate from the connective tissues and proliferation in the surrounding tissue to the difference of the decay rate and the rate at which fibroblasts are recruited to the thermal burn compartment. Thus, the steady state in the surrounding tissue for fibroblasts is represented by $[eqn]$

The steady state values calculated for the specific cell populations are then split between the undamaged and damaged variables dependent on the ionizing radiation dose. The proportion that is represented by the undamaged model variable is determined from linear quadratic cell survivability curves fit to experimental data [58–61]. The linear quadratic model was fit to experimental data using the fminsearch function in MATLAB using an objective function of the residual sum of squares between the model and the data and the default algorithm parameters for convergence. Calculations for the resulting initial conditions are show in Table 2 using the optimized parameters. The variable $[eqn]$ represents the radiation dose in Gy.

Table 2: Initial conditions for the variables in the surrounding tissue compartment.

The undamaged cell populations are initialized using the proportion of surviving cells and the damaged cells are initialized using one minus the proportion of surviving cells, as calculated by the respective linear quadratic model. State variables are assumed to be at steady state prior to and following irradiated thermal injury. Steady state values for macrophages ( $[eqn]$ ), lymphocytes ( $[eqn]$ ), and fibroblasts ( $[eqn]$ ) are strictly positive constants such that each is equal to its corresponding healed steady-state value.

Thermal burn compartment.

The thermal burn compartment includes the variables for damage and debris, the various inflammatory cells, fibroblasts, and pathogen. The damage and debris variables describe the severity of thermal injury and are initialized using thermal fluence alone. There is likely some radiation injury that occurs, but we assume that the thermal injury to cells causes complete cell death and thus supersedes any complete or partial DNA damage from ionizing radiation exposure. The population of inflammatory cells and fibroblasts are initialized to zero since we assume no viable cells will remain in the burn site.

Burn injuries considered in this model are caused by the direct absorption of radiant energy into the skin, referred to as flash burns [62]. The severity of the flash burn is used to initialize the damage and debris variables by linking the thermal fluence to a corresponding burn severity. To create a mapping of thermal fluence to burn severity, we estimated the deposition of radiant energy into the different layers of skin tissue, obtained a solution to the resultant bioheat transfer equation [10], and subsequently employed the Arrhenius damage equation [11] to determine the depth of severe cellular damage. Further details on this can be found in S5 File.

In our model, we associated the initial condition $[eqn]$ with a superficial thermal burn and $[eqn]$ with a superficial partial-thickness thermal burn. Due to the different cellular compositions of the epidermis and dermis, we developed a piecewise equation to model the level of tissue damage within the two layers. Equation 27 directly maps thermal fluence (f) to the expected damage constant for the equations for cellular damage and debris.

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and $[eqn]$ .

Pathogen is initialized using a combination of burn severity and radiation exposure. Radiation injury is associated with a risk of bacterial infiltration leading to infection [7,52,55] and we assume that more severe burns will provide greater opportunity for pathogen exposure from external sources. For ionizing radiation associated with concomitant wound infection, the organisms that cause most of the injury infections are skin flora organisms (staphylococci and streptococci) or environmentally acquired bacteria [63]. These bacteria are not significantly affected by ionizing radiation at the level considered in this model [64,65], so we assume that the ambient population is not reduced from the radiation environment.

In previous mathematical models featuring a pathogen population, the pathogen carrying capacity was typically used as the initial condition [29]. Risk and severity of infection-causing bacteria, however, varies directly with burn severity and radiation exposure [66,67]. We assume that there is no pathogen population when radiation exposure is not present, since superficial thermal burns alone are not associated with a significant risk of infection [66]. When radiation exposure does occur, we initialize the pathogen variable on a scale between zero and the carrying capacity where 0 is associated with no thermal injury and the carrying capacity, $[eqn]$ , is associated with the most severe thermal injury considered here, associated with a thermal fluence of 19 J/cm^2^. Thus, we define the initial conditions for $[eqn]$ using Equation 27 where $[eqn]$ is the initial condition for the pathogen variable, $[eqn]$ is the initial condition for the damage variable, and $[eqn]$ is the initial condition for the damage variable associated with a thermal fluence of 19 J/cm^2^. This scheme ensures that more severe burns are associated with a higher initial bacterial load and a prolonged pathogen presence. This acts as additional stress on the inflammatory system and corresponds with the increased immunesuppression and bacterial risk as burn severity increases [68].

$$[eqn]$$

Sensitivity analysis and parameter selection

The model contains a large number of unknown parameters which must be carefully selected to ensure biological feasibility. Due to the nature of the injury profile, data either does not exist or is extremely limited. Sensitivity analyses and leveraging of minimal data and biological trends were combined to select plausible parameter sets for the model. For sensitivity analysis, we employed the global method Multi-test Extended Fourier Amplitude Sensitivity Test (MeFAST) [69], which uses the Extended Fourier Amplitude Sensitivity Test (eFAST) algorithm [70,71] along with the t-test suggested in Marino et al. [72], the analysis of variance (ANOVA) with the Tukey procedure, and the Wilcoxon rank sum test utilized in Dela et al. [69].

Parameter combinations were initially generated using Latin hypercube sampling (LHS). For each parameter combination, the model was simulated for both a superficial thermal burn and a superficial partial-thickness thermal burn, as defined by the initial conditions of the state variables for damage and debris. Parameter combinations were assessed for three criteria, reflecting biological feasibility of the model predictions: (a) the simulation produced no numerical errors, (b) debris resolved before damage, and (c) the superficial thermal burn healed within $[eqn]$ days and the superficial partial-thickness burn healed within $[eqn]$ days [66,67,73,74]. The parameter minimum and maximum values for the eFAST algorithm were calculated from the sets satisfying these criteria.

To perform the sensitivity analysis, the model was simulated at five time points corresponding to days one, three, seven, ten, and fourteen post-injury. These are common data points used for collection in experimental models of thermal burn injury [75] and they correspond with identified response and peak times for neutrophils (1–7 days) and macrophages (3–14 days). Radiation levels were also varied using 1, 7, and 14 Gy. Sensitivity indices were calculated for varying values set for the number of samples per search curve ( $[eqn]$ ) and the resampling size ( $[eqn]$ ) to ensure the selection did not significantly influence the algorithm results. The largest values for $[eqn]$ and $[eqn]$ , 2000 and 32, respectively, were used in the final selection algorithm.

The sensitivity analysis results were used to reduce the estimated parameter space and aid in selection of a cohort of feasible parameter sets as well as a representative set of parameters that is indicative of the average model behavior. We developed an algorithm to select which parameters were considered highly influential (shown in S3 Fig). Those selected (shown in Table 3) were further sampled and outputs were assessed for biological feasibility. Non-influential parameters were set to their baseline value as determined by the average value of the initial sample range. Additional sampling details, including parameter inequalities that help to ensure stability of the healed state, can be found in S4 File.

Table 3: Parameters selected to be estimated or fixed using the sensitivity analysis results and the parameter selection algorithm.

Single parameter sets were solved numerically in approximately 0.05 seconds allowing us to simulate around 400,000 individual samples within a reasonable computational time and then assess the resulting transients for biological feasibility. A representative set was selected from the set of accepted sets by comparing the output of each parameter set to the temporal mean output of each model variable for the entire set of plausible parameter sets. The sets were compared quantitatively using both the weighted mean squared error between each set and the temporal mean as well as the coefficient of determination across all model variables.

Results

Model simulations

Model behavior was simulated using the representative set of parameters to demonstrate diverse model trajectories across initial conditions and outcomes associated with both the burn alone and the combined injury profile. The novelty of this model lies in the combined effects of burn with local irradiation, so various levels of radiation were explored. Simulations demonstrated delayed healing time, reduced immune cell infiltration, fibroblast dysfunction, and pathogen dynamics. All simulations use the initial conditions described in the methods section and include comparisons between the single injury (burn only) and the combined injury (burn and radiation).

Delayed healing times

An important feature of RCI is delayed healing for wounds [19–21]. This is captured by the model through increased healing times and prolonged presence of debris as shown in Fig 3. Both the 5 Gy and 14 Gy simulations exhibit increased duration of debris (Fig 3B) and slower resolution of damage (Fig 3A). The combined injury profiles also show a larger initial increase in damage as evidenced in Fig 3A by the initial positive slope. For the single injury profile, the initial positive slope resolves around 1.75 days, for the combined profile with 5 Gy, the slope remains positive until approximately 2.5 days and for injury with 14 Gy the slope is positive until around 3 days. This delayed healing is typically a result of delayed and prolonged immune cell infiltration at the wound site [22]. The prolonged presence of inflammatory cells also contributes to the initial increase in damage.

Plots of damage (Damtb) and debris (Debtb) for a superficial thermal burn with varying levels of radiation. Panel A shows Damtb and panel B shows Debtb. Superficial thermal burns alone should heal fully within 4 to 5 days. The introduction of ionizing radiation delays the immune response and results in prolonged healing times.

Reduced immune cell infiltration

Local radiation injury is known to delay the infiltration of immune cells from the bloodstream and into the wound due to the damage of microvasculature in the affected tissues [24]. This is included in the model by inhibiting influx rates into the surrounding tissue from the bloodstream which results in reduced resource pools available for activation at the wound. Fig 4 shows the cells in the surrounding tissue supplied by the populations in the bloodstream. As shown, combined radiation injury results in reduced initial levels of undamaged cells (due to the initial radiation injury) which is further inhibited by microvasculature damage in the early stages of inflammation (within the first few days). These delays propagate to the corresponding activated populations at the thermal injury site.

Plots of undamaged cells in the surrounding tissue supplied by the bloodstream.Panel A represents neutrophils, panel B represents macrophages, panel C represents lymphocytes, and panel D represents fibroblasts. Populations in the surrounding tissue are affected by both the initial cellular injury from ionizing radiation as well as reduced influx due to microvasculature damage in the surrounding tissue.

For the combined injury profiles, all immune cells in the surrounding tissue remain below the single injury levels, aside from neutrophils. This population returns to zero once the wound is healed; since healing is delayed in the combined injuries, levels remain elevated. An initial dip in the cell population in the surrounding tissue is expected, like those shown in Fig 4B and 4C, since immune cells will be attracted to the injury. In the combined injury profiles, extravasation from the blood stream is necessary to support both activation of immune cells in the wound as well as replacement of the damaged resting population.

Similar trends can be observed in the thermal burn compartment (Fig 5). All variables show marked delays in immune cell peaks as well as reduced magnitude of response when radiation is present. M1 macrophages and L1 lymphocytes actually exhibited some differing behavior with peaks occurring slightly earlier with radiation exposure of 14 Gy. This can be a result of the radiation exposure encouraging an increased inflammatory environment and delaying the switch to anti-inflammatory phenotypes represented as M2 macrophages and L2 lymphocytes. Both are significantly delayed with the 5 Gy and 14 Gy exposures.

Plots of immune cells in the injury compartment.Panel A represents neutrophils, panel B represents M1 macrophages, panel C represents L1 lymphocytes, panel D represents M2 macrophages, and panel E represents L2 lymphocytes. Immune cell response in the injury is delayed due to ionizing radiation exposure.

Fibroblast dysfunction

Ionizing radiation affects dermal fibroblasts by reducing proliferative abilities [23]. Fibroblast dysfunction, along with reduced infiltration as shown in Fig 6, results in reduced fibroblast populations and irregular collagen deposition. As shown, these effects increase with ionizing radiation levels and correlate with delayed healing times.

Plots of fibroblasts (Ftb) and damage (Damtb) in the thermal injury.Ionizing radiation exposure alters fibroblast function resulting in reduced proliferative abilities and irregular collagen deposition, delaying damage resolution.

Sustained pathogen population

RCI is associated with an increased risk of bacterial infiltration and sustained infection [2,7]. Superficial thermal burns, on the other hand, are not typically associated with a risk of infection [66]. As shown in Fig 7, the model has the capacity to induce prolonged pathogen presence that could result in a systemic infection if spread to the bloodstream. Without resolution of the pathogen variable, the thermal injury will fail to heal, resulting in chronic inflammation.

Plots of pathogen (Ptb) and damage (Dam_tb) at the thermal injury site.Since ionizing radiation is associated with risk of bacterial infiltration leading to infection, the model includes the capability to capture prolonged bacterial presence and delayed healing which could result in a non-healing wound and chronic inflammation.

Discussion

We developed a model for a combined injury profile consisting of a superficial thermal burn and locally confined ionizing radiation exposure. The model consists of a system of ODEs representing key inflammatory cells and fibroblast populations in the thermal injury and surrounding tissue compartments. A mechanistic model for combined injury provides time-dependent outputs that can capture the synergistic effects of combined injury for a wide range of inputs while also offering pathways for further predictive capabilities associated with the injury. This model serves as an initial exploration into modeling combined injury using a mechanistic model. A mechanistic model presents significant advantages over previously developed empirical models that are only able to predict single variables, such as injury prevalence [27] or probability of survival [28]. Additionally, mechanistic models can be used to better understand the underlying mechanisms involved in RCI, such as specific immune cell dynamics and pathogen dynamics. This model will be further developed to improve predictions of first responder capabilities in a nuclear incident.

The formulated model presents a mechanistic modeling approach, which relies on known biological mechanisms involved in inflammation and early proliferation, as well as experimental observations and assumptions from animal trials and first responder reports from historical events. Due to the nature of the injury, real world data is not available. Despite this, the model still uses a number of qualitative assessments to ensure biological feasibility of the results. The simulations of the single injury profile (burn only) were restricted to those that healed within the previously determined healing times for superficial and superficial partial thickness thermal burns [66,67,73,74]. The model formulation and simulations also rely on known dynamics for inflammation and activation, including both timing of peaks in specific cell populations and interactions between the cells. Observed behavior in the combined injury profile was also controlled for in parameter selection (S4 File) and verified qualitatively in model outputs. While this does not provide quantitative measures of model verification, this form of modeling can still offer useful predictions about the biological behavior that can be further improved if data becomes available for the model outputs or model parameters. The present model serves as the first iteration of this type of modeling which will continue to be improved.

The model was able to recreate various known dynamics of RCI with burn including increased healing times, reduced immune cell infiltration, fibroblast dysfunction, and sustained pathogen. Experimental models for RCI, both with burn [8,76,77], and other injuries (laceration) [78–80], consistently report prolonged healing times and reduced immune cell activity. Fibroblast dysfunction, including both reduced proliferation and irregular collagen deposition, has also been evidenced in the experimental studies [23,51] and further contributes to delayed and deficient wound healing.

Pathogen dynamics are particularly important in combined injury models involving burn injury. In superficial thermal burns alone, there is little to no risk of infection [66,67], but ionizing radiation is associated with a risk of infection regardless of whether there is a localized [52] or whole-body exposure [53–56]. This risk is further elevated for irradiated burns [7]. The majority of accepted parameter sets exhibited quick resolution of any foreign infiltrates, as would be expected for less severe burns, but the model also contains the capacity for a sustained pathogen response (Fig 7). The pathogen equation (Eq. 25) accounts for minor immunosuppression during radiation exposure from reduced availability of macrophages due to cellular death. This, however, has not been accounted for in the background immune response, which is mainly mediated by mast cells and natural killer cells [81]. For this reason, the current formulation may not accurately predict sepsis risk and should be revisited in future iterations and when considering burns with larger total body surface area and/or penetration depth or those with whole-body radiation exposure. Since infection observed in radiation injury could also be the result of bacteria entering the blood stream from damage to the gastrointestinal tract [82], it is important to continue to consider pathogen dynamics when models of whole-body radiation exposures are considered for RCI.

Our modeling approach does have some limitations. First, this initial model formulation contains some necessary assumptions to support modeling, including instantaneous exposure to thermal and ionizing radiation at the injury site and some cell based assumptions. The present model serves as the first iteration of this type of approach for combined injury and will continue to be improved upon. Second, the model is confined to the local domain. This prevents the model from representing more severe injuries that result in systemic effects [3,76] and whole-body radiation which results in hematopoietic disruption and additional systemic effects [18,82]. Third, the model does not account for additional injuries, such as lacerations or fractures, which can put additional stress on the body, potentially exacerbating immune response delays. Fourth, the model does not account for specifics of the thermal burn including location and percent total body surface area. Since we assume the thermal injury does not initiate a systemic response, percent total body surface area is assumed to be relatively low. Location is also important in determining burn severity since burns are typically more severe on the face and hands.

Despite its limitations, the model can still offer useful applications. One application we’ve explored is using immune cell levels at the wound site to predict medical outcomes, namely hyperalgesia. The model variables describing the inflammatory response can be used to predict the experience of inflammatory pain by estimating the biological concentration of signals detected by nociceptors, sensory neurons that respond to damaging stimuli. Inflammatory pain is associated with tissue damage where signaling biomolecules are released from damaged cells and inflammatory cells which prime nociceptors and increase the sensitivity of the nerve endings [83]. The signaling molecules causing an inflammatory pain response and hyperalgesia include IL-1β, TNF-α, and IL-6, which are mediated and released by active inflammatory cells. Model variables can estimate this response and estimate the occurrence of pain in RCI profiles.

To further develop the combined injury framework, future models should extend to systemic injuries including both more severe burn injuries and widespread radiation exposure. When considering widespread radiation, the hematopoietic system and gastrointestinal tract are particularly sensitive and should be included in whole-body exposure models. Systemic injuries would also allow for the estimation of systemic symptoms affecting medical outcomes. This could include fatiguability and weakness which largely result from increased circulating cytokines [84] and upper-gastrointestinal distress which results from major cell death in intestinal crypts [82] and serotonin release in the gut [85–87]. Mechanistic models for RCI can be further developed to include these systemic effects and serve as sophisticated modeling tools to estimate casualties and capabilities of first responders in the event of nuclear incident.

Conclusions

Radiation combined injuries present a formidable modeling problem since single injury and empirical models are unable to properly capture the synergistic effects associated with concurrent injuries and provide accurate extrapolatory predictions. We developed a mechanistic model of RCI including local radiation exposure and burn injury which allows for estimation of injury effects outside the restriction of available data. A mechanistic model also provides a means to explore medical outcomes associated with the inflammatory response that could affect first responder performance. This model serves as a basis for developing mechanistic combined injury models of systemic injury associated with whole-body radiation exposure and more severe injuries such as larger burns and major fractures. Further development to include systemic effects would improve estimations, providing tools for medical planning and triage in a nuclear incident.

Supporting information

S1 TableParameters with descriptions and selected values.Model parameters including descriptions and units. The listed numerical value was determined by the representative set.(PDF)

S2 FileDetailed description of model equations.(PDF)

S3 FigFlowchart of parameter selection algorithm to determine whether a parameter should be fixed or further sampled.(TIF)

S4 FileParameter sampling and acceptance criteria.(PDF)

S5 FileMapping thermal fluence to damage and debris.(PDF)

S6 FileSensitivity analysis results.(PDF)