The Modelling of the Post-Operative Perfusion in Burns from LDI Data

The paper deals with an empirical mathematical model serving for the time evolution of postoperative blood perfusion in burn wounds developed in the setting of surgical treatment of full thickness, or non-healing deep dermal burns involving application of Split Thickness Skin Grafts (STSG), and Autologous Platelet Concentrate (APC), with the main goal of early assessment enhancement and healing prognosis of burns. Numerical parameters of the model including a 95 % confidence level prediction interval were determined using the nonlinear fit procedure. We cautiously claim that the suggested model could serve as a beneficial tool in objective comparisons of treatment efficacy and facilitating the decision making process. However, reliably determined model parameters still require a large amount of laser Doppler imaging data to be processed for individual kinds of treatment.


Introduction
Although Laser Doppler Imaging (LDI) has been entrenched and almost ubiquitous in burns diagnostics since it first became available, one can observe persistent striving for new means of LDI data processing which would be able to contribute additional benefits, improve the result accuracy and reliability, set forward objective comparison of treatment efficiency [1], [2], [3], [4], [5] and facilitate local treatment strategy decision making [6], [7] and [8].The aforementioned reasons, together with increasing emphasis placed on outcomes following treatment of deep dermal and particularly full thickness burns, which are undergone to surgical treatment, delivered the main impetus for the development of a novel tool and inspired the idea of its advantageous embodiment into a simple yet powerful mathematical model of post-operative perfusion time evolution with the aim of refining both instantaneous burn assessment and burn healing prognosis [9], [10] and [11].
As a background and source of data for setting up the mathematical model, an ongoing research of surgical intervention involving application of autologous platelet concentrate and Split Thickness Skin Grafts (STSG) was adopted.Autologous Platelet Concentrate (APC) is plasma with a concentration of blood platelets elevated above baseline [10], [12], [13] and [14].

Methods
The study protocol was reviewed and approved by the Institutional Review Board of the University Hospital Ostrava.All patients were fully informed about the study, and gave informed consent to the treatment and the subsequent measurements.
The study presented is a non-randomized prospective monocentric data analysis based on objective testing of outcomes of surgery on deep burns.All procedures involved in the study were carried out by the objective measurement of blood perfusion using the laser Doppler imaging technique (see Ref. [15] and references therein).

Patient Criteria
From March, 2009, to November, 2012, a total of 27 patients (15 men and 12 women) underwent surgical treatment of full thickness and non-healing deep dermal burns with application of STSG and APC.
The criteria for including the particular patient were as follows: • an adult patient (at least 18 years of age) with full thickness burns and/or non-healing deep dermal burns, • a patient who signed the informed consent form prior to participation, • absence of hematological malignancies, chronic renal failure, inborn thrombocytopenia, current intake of oral anticoagulants, platelet aggregation inhibitors or abnormal bleeding history.

Material and Instrumentation
The APC and autologous thrombin were prepared by density gradient centrifugation of patient's peripheral venous blood collected using standard venipuncture techniques, most frequently from an upper extremity peripheral vein using the Harvest SmartPReP Platelet Concentrate System (Harvest Technologies Corporation, Plymouth, MA, USA).The blood was collected before the commencement of surgery, the APC preparation and necrectomy of the thermally destroyed skin proceeded simultaneously.The blood was centrifuged in accordance with the manufacturer's recommendations.
The necrotic area was then grafted with STSG.Upon application of APC on STSG-transplanted areas and activation of thrombocytes by autologous thrombin, the platelet α-granules granulated releasing growth factors acting as healing stimulants [16].The entire process was carried out under strict sterile conditions in the operating room, simultaneously with the surgery.
For evaluation of blood perfusion in early phases (less than 14 days after the surgery), the PeriScan PIM 3 LDI scanner by Perimed AB, taking advantage of the principles of laser Doppler imaging, was used [17], [18], [19] and [20].The system was accompanied with respective software LDPIwin 3.1 supplied along with the hardware.Figure 1 shows the freehand thick white closed curves defining the RoIs were drawn with regard to a reasonable representation of the burn areas yet to keep clear of the inflammatory periphery which would introduce falsely elevated perfusion values.The mean perfusion within the area encircled by the curve was figured out using the LDPIwin 3.1 program.(These grayscale images were converted from the LDPIwin 3.1 program's original color-coded images using a utility written in the system Mathematica.

Perfusion Data Acquisition
Laser Doppler scanning was carried out during each dressing change, first before the surgery, subsequently on every even Post-Operative Day (POD), starting on the 2 nd and ending on 14 th day.
All scans were conducted according to the same preset protocol: on each patient at the same time of day, in the same examination room with controlled room temperature between 22-24 • C and relative humidity of 50 % ± 10 %.The ambient light conditions were ensured to comply with the recommendation of the scanner manufacturer [21], and the device was calibrated regularly using the calibrating phantoms.All patients were inactive for a period of 15 minutes, and were always examined in the same position.The patients were asked to remain at rest to prevent movement artifacts.The scanning head was adjusted to a distance of 15-20 cm (about 6-8 in) from the examined area, the axis of the measuring head was being set as perpendicular to the skin surface as possible.The scanner settings were recorded and saved for future reference.
The resulting 2D perfusion maps were preprocessed using the LDPIwin 3.1 program as illustrated in Fig. 1.All operators introduce idiosyncratic errors, in order to avoid varying errors caused by the operator.
Due to the demands of clinical situations, only 9 patients out of a total 27 passed through all seven scheduled measurements.Consequently, some datasets are incomplete the only complete ones are from PODs 6, 8 and 14.

Convectional Statistical Analysis
Seven perfusion datasets from consecutive postoperative days 2, 4,...,14 were primarily undergone to conventional statistical analysis with the particular intention to lift the veil off its pitfalls.All statistics (including the modeling approach calculations to be discussed later) were carried out using the system Mathematica [22] and [23], a significance level of 5 % was adopted.
The datasets were tested for normality and equal variances.The null hypotheses that the datasets were drawn from a normally distributed population, and that the datasets had equal variances, were found unlikely and could not be retained.Thus, the datasets failed to meet basic pretest assumptions of parametric multi-comparison testing [24].
We run into problems, however, with independence because the datasets are generally of different sizes.Three approaches to tackle the problem could be adopted: • taking into account completely measured patients only, • using linear regression substitution for missing perfusion data, • performing post-hoc regression substitution for missing perfusion data using the model developed for and applied to individual cases.
Neither of approaches 2 nor 3 is exactly "missing data imputation" in the sense used by, e.g., Ref. [9], but for brevity we retain the term hereinafter.While the number of complete datasets drops down to 3 (with the number of patients measured in all 7 days being 9), using imputation enables the full number of datasets to enter the independence tests at the cost of imputationintroduced bias.
Testing under the first approach retained the null hypothesis that no association exists for any of 21 dataset pairs.Testing using the same tests under the second approach led to independence refusal by at least three tests in nine pairs of POD.The situation did not change much when testing under the third approach, although it takes advantage of a more exact regression model, independence was refused by at least three tests for seven pairs.Taking the above findings into account, a base 10 logarithmic transformation was applied to approximate the normal distribution more precisely.Transformation is also a good remedy for heteroscedastic datasets and their grossly positive skew.
Repeated testing for normality revealed that p-values of all tests in all post-operative days do not drop below 0.11 (with typical value above 0.4), the distributions of original data can be deemed log-normal.Testing for equal variances does not reject the null hypothesis that the datasets have equal variances.In the following, transformed data will exclusively be used without explicit emphasis.
Figure 2 shows graphically expressed results of the conventional statistical analysis of transformed perfusion data.It can easily be observed that both perfusion medians and means monotonically decrease through 2 nd to 8 th post-operative days, reaching stagnation value in the remaining days, perhaps with the exception of POD 14.The most expressed dynamics is noticeable in early post-operative stages, as might be expected.
Based on pretests for normality and homoscedasticity results, the transformed perfusion data are qualified for testing for equal means or medians in individual postoperative days.A single factor analysis of variance (ANOVA) was performed on the data resulting in p = 0.00025.Subsequent post hoc tests indicated significant difference in the following pairs of post-operative days: (2, 8), (2,10), (2,12), (2,14), and (4,14).
Because of problematic independence of datasets, imputed data were tested using Friedman rank test for paired data.For both linear and post-hoc model regression imputation, the results induced rejection of means equality (the p-value was of order 10 −10 or less).Subsequent post-hoc testing indicating significant differences in the following pairs in addition to those revealed by ANOVA post hoc tests: (4, 6), (4, 8), (4, 10), (4,12), and (4, 16).

Building the Model
Let us remind that, as shown in Fig. 2, the perfusion means and medians monotonically decrease from 2 nd to 8 th post-operative days, reaching stagnation value after that.Consequently, in order to capture the nature of the perfusion behavior, we can exclude polynomials of arbitrary nonzero order (including of linear) as the model function, since they diverge to plus/minus infinity for large time.The simplest and most natural way of portraying the perfusion drop features is provided by the exponential formula: where t denotes time (in units of days) elapsed since surgery, P s stands for the steady perfusion value for large time t, P 0 is the difference between (hypothetical) elevated perfusion value at t = 0 and the steady state (in other words, P s + P 0 is the perfusion value at t = 0), and λ is another parameter expressing the rate of restoration toward the steady state.Exponential formulas like Eq. ( 1) are characteristic of a broad class of natural processes which are gradually steered to a steady, unchanging state.Qualitative behavior of such model is shown by the leftmost panel in Fig. 3. Denoting the difference between the immediate perfusion and the steady state perfusion value, and the rate of return to steady state, respectively, a simple property of Eq. ( 1) implies Ṗ = −λD, i.e., the immediate rate of return is proportional to the immediate excess of perfusion above the steady value.The absolute value of the rate of return P reaches its maximum at the very beginning, t = 0, and converges to zero for late post-operative days, Ṗ (0 This behavior is illustrated in the leftmost panel of Fig. 3 by the downward tilted dashed tangent.The constant λ can be referred to as "recovery constant" related to a more suitable characteristics of the decaythe recovery halftime -by the formula T 1/2 = 0.693/λ.The meaning of the recovery halftime is that its every elapse reduces the difference D in half.Of course, the need not be taken seriously for time from t = 0 (surgery day, POD 0) to t = 2 (POD 2).
Nonlinear regression analysis using the model Eq. ( 1) applied on the transformed data indicated that determination of the key coefficient λ is statistically insignificant, yielding p-value above 0.05 for its estimation, or, stated otherwise, its interval of confidence contains zero (see the rightmost panel in Fig. 4).This leads to a proposal to modify the model Eq. ( 1) by introducing an exponent γ > 1 as follows, Figure 3 shows reasoning for the model selection.Purely exponential function P (t) describing the restorative process displayed in the leftmost panel shows the steepest gradient at the very beginning (t = 0), i.e. the tangent (dashed arrow) has γ the most negative slope at that point.For any γ > 1 the function P (t) brings an essential change in its behavior, turning the initial tangent to the horizontal direction.With the parameter γ rising (while keeping the other parameters constant), an increasingly apparent plateau, delaying the decline to a later phase, is progressively formed.This is however compensated by expediting the subsequent decrease and by earlier establishing a steady state.These features (which bestowed the name "suspended/expedited Gaussian recovery" upon the model) are qualitatively demonstrated in the remaining three panels for γ = 4/3, 3/2, 2 and 2.5, respectively.The shaded area delimited by the curve, perfusion ordinate, and line at the steady state level parallel to the abscissa, is proportional to the excess of total number of blood elements pervaded into the wound per unit area during entire healing period.
Figure 4 shows incrementing the parameter γ by 0.1 from 1.0 to 3.0 (and adding some extra values of small integers ratio like 6/5, 5/4, 4/3, 5/3, 7/4), a nonlinear regression using Levenberg-Marquardt method was carried out for fitting model parameters on transformed perfusion data.The panels show the dependence of the p-values for parameters P s , P 0 and λ (from left to right) on γ.While the dependence of p-value of P s and P 0 is monotonically decreasing, the p-value for parameter (λ) minimizes at (γ) equal to 2.0 or very close to 2.0.But this means that at that value the 95 % confidence interval minimizes its width.Best fit in the sense of minimized p-value of the most critical parameter λ was achieved for γ = 2, as shown in the rightmost panel.Figure 5 shows the points representing the deviations (residuals) of transformed perfusion values from those predicted by the best fit parameters of model Eq. ( 4).The residuals are evenly distributed for all values and show no perceptible trend.The qualitative properties of this model are similar to those of Eq. ( 1) with the exception that the initial drop rate vanishes, Ṗ (0) = Ṗ (∞) = 0.
The effect of the new parameter γ > 1 is demonstrated in Fig. 3.While for γ = 1 the tangent at the initial point is pointing askew, for γ > 1 the initial tangent is always horizontal, developing small plateau in the vicinity of the origin.The larger γ is, the more the plateau extends to the right.The delay in decrease is then compensated by expedited drop and quicker settling on a flat steady state level.
Let us now dwell on developing the model based on Eq. ( 3).The nonlinear Levenberg--Marquardt regression method was used for fitting the model parameters on transformed perfusion data.Technically, the regression makes use of the fact that the model log 10 P (t) = log 10 (P s + P 0 e −λt γ ) takes advantage of favorable properties of transformed data.The fitted parameters were P s , P 0 and λ, while an ansatz was used for the exponent γ incremented by step of 0.1 from 1.0 until 3.0.The best fit, the results of which are summarized in Tab. 1 and Fig. 3, was achieved for γ = 2 (details are explained in figure caption) with adjusted coefficient of determination R 2 = 0.986.
Tab. 1: Results of regression analysis on the model Eq. ( 1).
The 68 % (95 %) confidence interval corresponds to ≈ 1σ (2σ) departure from the best-fit estimate.The adjusted coefficient of determination is R 2 = 0.986.Finally, the results were verified by means of analysis of residuals.The residuals in Fig. 5 show no perceptible trend.Testing for normality using Anderson-Darling, Shapiro-Wilk and Jarque-Bera ALM tests yields p-values 0.34, 0.43 and 0.77, respectively, so the normality of residuals has to be retained.Moreover, the Student t-test effectively yields p = 1 (sic!, the value of the t-statistic −2.53 • 10 −15 ) for the hypothesis that the residuals mean equals to zero.Thus, the model Eq. ( 3) with γ = 2, appears to be highly compliant to the data under investigation.
There are, however, some points worth clarifying.The model function Eq. ( 4), rewritten here for convenience with numerical best-fit results substituted from Tab. 1 seems to provide reasonable foundation for processing of the perfusion data collected in the way described.Looking at the graphic representation of formula Eq. ( 5) in Fig. 6(a) (confidence level 95 %) and Fig. 6(b) (confidence level 68 % corresponding to standard deviation ≈ 1σ), the perfusion drop is clearly manifested up to roughly 10 th post-operative day.Starting from POD 12 the steady state 91 PU is virtually reached.For comparison, the results for POD 30 are marked with a circled dot -it is located at the borderline of the 95 % dark gray model reliability band, but it has nothing to do with the prediction of a single perfusion value, which is marked by the light gray prediction band.
There are some points that distinguish the model described by Eq. ( 4) from that described by Eq. ( 1).While the latter permits the perfusion to drop from the very initial time (immediately after the surgery), the former delays the drop for a while at the beginning, as if the perfusion drop was inhibited by an unknown factor that later vanishes, replacing delay by hastened decrease (cf.Fig. 3).It should be admitted that the nature of that lag remains undisclosed.One need not to take a grave view of the model prediction for time between t = 0 (surgery day, POD 0) and t = 2 (POD 2).On the other hand, within the interval of monitoring (starting from the first post-surgery month) it is capable to reflect the real progress of healing process reasonably.The early perfusion plateau predicted by the model can be an artifact resulting from the fact that the model fails in description of the perfusion evolution in early stages after the surgery (although it may be fairly good approximation from POD 2 to 14), or it can harbor rudiment of still concealed phenomenon.However, it is too early for such a complex issue to be judged -but provided Eq. ( 4) turns out to model the perfusion reasonably from POD 2 on, the early behavior surely deserves further investigation.
Interestingly, there is another theoretical aspect to be emphasized.The rate of return to steady state Eq. ( 2) for model dependence Eq. ( 4) reads Ṗ = −2λtD.
Figure 6 shows the result of perfusion time evolution nonlinear fit.The perfusion axis uses logarithmic scale.Model Eq. ( 4) was used with parameters P s , P 0 and λ (the meaning of individual parameter is explained in text).Large black points represent mean values of transformed perfusion in post-operative days, small gray points represent individual transformed perfusion values (horizontal random dispersion is deliberately introduced to prevent overlapping of close values).The 95 % prediction band of reliability for future measurements is in light gray, white curve in the middle represents the best fit, with a dark gray 95 % model reliability band alongside.The vertical lines indicate the critical recovery time tc and its uncertainty.The circled dot shows the position of the perfusion mean on the 30 th post-operative day. Figure 6(b) same as in Fig. 6(a) but with 68 % prediction band of reliability for future measurements (light gray) and model reliability band (dark gray).The 95 % confidence level is however retained for uncertainty of critical recovery time t c , its narrowing resulting from change of reliability band margins.(The value 68 % corresponds to ≈ 1σ interval of confidence.)Although the model formula Eq. ( 4) was -and its parameter γ = 2 in particular derived purely on empirical basis, the recent formula offers a simple and plausible interpretation: the rate of decrease is proportional both the difference between immediate and steady state perfusion and the time elapsed from surgery.On the other hand, this reasoning lacks deeper insight and theoretical background.We present it here mainly for completeness and as an initiative for further research.
The conformity of the model is at the expense of the concept of recovery halftime indeed, the parameter γ > 1 causes gradual decrease in subsequent "halftimes."However, remembering the nature of perfusion, even better estimation of recovery time can be introduced: the left-hand side denominator of expresses the excess area between the perfusion curve and the steady state horizontal exhausting entire light gray plus dark gray area under the curve in Fig. 6(a) is a measure of total blood perfusion during healing per unit area, while the numerator quantifies the same but only up to certain critical "recovery time" t c .Comparing their ratio, as shown in Eq. ( 6), to a suitable chosen cutoff value L c , typically L c = 0.95, after some manipulations using Eq. ( 2) and Eq. ( 4), we arrive to an explicit expression for the recovery time t c = erf −1 (L c ) √ λ , where erf −1 denotes inverse error function tabulated in all standard references like [14] and [25].For L c = 0.95 this formula simplifies to t c = 1.386 √ λ .
The recovery time can be interpreted as the period during which 95 % of total excess perfusion takes place.Substituting the estimate of the parameter λ provides t c = 8.0 days, which is value consistent with results obtained using conventional statistical analysis.The t c 's confidence interval was calculated numerically from the marginal model curves; the results are shown in Fig. 6(a) and Fig. 6(b) for parameter confidence level 95 % and 68 %, respectively.

Conclusion
We cautiously claim that the suggested model could prove itself to be a beneficial tool in efficacy comparison of treatments under research as well as in fast and reliable routine clinical evaluation of post-operative burn condition, thus facilitating decision making.However, reliably determined model parameters still require creation of a large LDI database to be processed separately for various kinds of treatment, with the main aim of establishing reference bands for reliable detection of perfusion behavior, its probable extrapolation, and subsequent decision.The early detection of deviating perfusion (and/or its time change) implies a possibility of urgent healing measures in sufficiently early phases when an increase of complications has not yet been clinically manifested.The presented model predicts perfusion plateau in the very early stages after surgery, however, future studies are needed for deeper insight into its origin and biological role in burn healing.

Fig. 1 :
Fig. 1: An example of a region of interest (RoI) selection.

Fig. 2 :
Fig. 2: The box-and-whisker chart for logarithmically transformed perfusion in seven PODs examined.The thin white lines indicate medians, short thick ticks indicate means and circles mark outliers.

4 :
Modeling of the p-values for parameters Ps, P 0 and λ (from left to right) on γ.

Fig. 6 :
Fig. 6: The result of perfusion time evolution nonlinear fit.