Optimal switching time control of the hyperbaric oxygen therapy for a chronic wound

1.   Introduction

Chronic wounds, which fail to heal through the normally orderly process of stages and remain in a chronic inflammatory state, are a major burden for health care systems and a significant socioeconomic problem worldwide ^[1,2]. For example, 2.4–4.5 million people are believed to suffer from the chronic lower extremity ulcers (a typical chronic wound disease) per year in the United States ^[3]. These ulcers last on average 12 months and recur up to $60\%$ to $70\%$ of patients. In a Dutch nursing-home population, a recent study shows a $4.2\%$ prevalence rate of chronic wounds ^[4]. Chronic wounds lead to loss of basic functions of human body and decrease quality of life tremendously. Along with other diseases such as diabetes, $24\%$ of these wounds will ultimately result in amputation of the lower limb. Thus, effective treatment for chronic wounds is urgently needed.

The dynamic evolution process of the chronic wound mainly consists the oxygen, neutrophils, invasive bacteria, and chemoattractant (a chemical that attracts neutrophils to the bacteria). There are many useful mathematical models that address these four physical and chemical phenomenon. McDougall et al. model the dermal wound healing by a weighted average of cells' previous motion direction ^[5]. Cumming et al. use a cell-based discrete approach to describe the delicate interaction between elements that vary widely in nature and size scales ^[6]. Gaffney et al. use standard form of the conservation equation for the tip concentration and the endothelial cell density ^[7]. Thackham et al. consider advection-dominated wound healing models with a source term and obtain an analytic solution under some simplification ^[8]. Vermolen has developed a system of two dimensional nonlinear reaction-diffusion equations, in which oxygen, growth factors, epidermal cells are considered ^[9]. Essentially, the dynamic evolution process of the chronic wound is spatiotemporal evolution, which is generally modeled by coupled partial differential equations (PDEs). For example, continuum reaction-diffusion models, one of typical PDEs, are the most commonly theoretical approach to study the angiogenesis process ^[10].

Wound healing is a highly regulated and undeniably complex process, which can be divided into four interconnected and overlapping stages: haemostasis, inflammation, proliferation and remodelling ^[11,12]. Among these processes, hyperbaric oxygen therapy is a critical means for healing the wound and has been shown to improve the healing of chronic wounds in clinical trials ^[13]. Hyperbaric oxygen therapy is breathing 100% oxygen under increased atmospheric pressure. The primary mechanism behind the use of this therapy for wound healing is to correct the amount of oxygen delivered to the wound site. The positive effects stem from increased tissue oxygen tension and pressure within the wound site, such as increased collagen synthesis, modulation of the immune system response, inhibition of bacterial propagation, prevention of leukocyte activation and so on ^[14]. Kalliainen et al. conclude that topical oxygen has no bad effects on chronic wounds and show the beneficial indications in promoting wound healing ^[15]. Gajendrareddy et al. demonstrate that systemic oxygen therapy ameliorates stress-impaired dermal wound healing ^[16]. Eisenbud gives the correction of wound hypoxia and clarifies that oxygen is beneficial to many aspects of healing but it does not necessarily follow that the more is better ^[17]. These work described above greatly promotes the hyperbaric oxygen therapy and provides guidance for the medical staff. However, these results are from a large number of repeated experiments and feel weak for individual difference treatment. Moreover, Gill et al. discuss the application, mechanisms and outcomes of the hyperbaric oxygen therapy but also mention that the therapy is expensive ^[18]. Thus, this paper considers the hyperbaric oxygen therapy for a chronic wound from the perspective of control synthesis and economic cost.

To our best knowledge, there are few literature studies in the control of the hyperbaric oxygen therapy for a chronic wound. The time and amount of oxygenation applied during the hyperbaric oxygen therapy provide the regulatory mechanism for the advanced control strategy study to improve wound healing ^[19,20]. Recently, Guffey demonstrates the existence of the state solutions for a bacterial infection model with modified boundary conditions ^[21] and establishes some theoretical foundations of optimal control in a chronic wound. However, the convergence of the state system forward and the adjoint system backward is not solved and the implementation of optimal control is blocked. In this paper, we complement the optimal switching time control of the hyperbaric oxygen therapy for a chronic wound based on the sensitivity analysis and have a great prospect in the application and guidance of the hyperbaric oxygen therapy.

The paper is organized as follows. In Section 2, the dimensionless oxygen concentration, dimensionless neutrophil concentration, dimensionless bacterial concentration and dimensionless chemoattractant concentration with corresponding boundary conditions and initial conditions in a chronic wound are modeled. Then, the optimal control problem for inhibiting the growth of bacterial concentration through the hyperbaric oxygen therapy is proposed. In Section 3, The method of lines, which only discretizes the space while keeps the time continuous, is used to reduce the PDE optimization problem into an ODE optimization problem. In Section 4, the control design problem of the hyperbaric oxygen therapy is the parameter selection problem of the optimal time switching points essentially. However, the switching time problem is difficult to solve. Thus, the time-scaling transformation approach which turns the original problem into a new problem with fixed switching time is applied and the gradient of the cost functional corresponding to the time intervals is derived based on the sensitivity analysis. In Section 5, computational numerical analysis is carried out to show the effectiveness of the optimal switching time control for bacterial concentration suppression under different discretization numbers of the time intervals. Finally, Section 6 concludes the paper with comments and suggestions for future research.

2.   Problem formulation

We consider the situation, where the chronic wound is considered to be one-dimensional with the spatial distance $x = 0$ located at the center of the radial wound site and $x = 1$ located at the edge of the chronic wound nearest healthy dermis. This paper mainly uses oxygen, neutrophils, bacteria and chemoattractants to describe the dynamic evolution process of a chronic wound. Moreover, the relationship of these variables' interaction is shown in Figure 1.

2.1.   Dimensionless oxygen concentration model

The movement process of the dimensionless oxygen concentration $O(x, t)$ in a chronic wound is governed by the following PDE model ^[22,23]

where $t$ is the time variable; $\eta$ is a constant input of external oxygen; $\kappa$ is the oxygen coefficient of the hyperbaric oxygen therapy; $\lambda_{NO}$ is per unit neutrophil; $\lambda_{BO}$ is per unit bacteria; $N(x, t)$ is the dimensionless neutrophil concentration; $B(x, t)$ is the dimensionless bacterial concentration; $\lambda_{O}$ is the natural decay rate of oxygen. Oxygen lost in any other way, unrelated to the neutrophils and bacteria, occurs at a constant rate of $\lambda_{O}$. For example, other leukocytes need to consume oxygen. Oxygen naturally decay is also described in ^[24]. Moreover, during the hyperbaric oxygen therapy, we can only choose the oxygen tank on or off. Thus, $u(t)$, a time-dependent oxygen therapy, is given as

Since the oxygen concentration satisfies the natural boundary condition at the center of the chronic wound and is enough at the edge of the chronic wound nearest healthy dermis, the boundary conditions for (2.1) are:

We assume the initial state of the oxygen concentration is enough and then obtain the following initial condition for (2.1)

2.2.   Dimensionless neutrophil concentration model

The movement process of the dimensionless neutrophil concentration in a chronic wound is governed by ^[25,26]

where $D_N$ is the diffusion coefficient of neutrophils; $\kappa_N$ is the attraction coefficient of neutrophils via chemoattractant gradient released by the bacteria; $C(x, t)$ is the dimensionless chemoattractant concentration; $\alpha_{BN}$ is the proliferation coefficient of neutrophils corresponding to the bacteria and oxygen; $\lambda_N$ is the decay coefficient of neutrophils dependent on the presence of bacteria and oxygen, $e$ is the relative coefficient of neutrophils corresponding to the bacteria. As bacteria proliferate in the wound, neutrophils are recruited into the chronic wound site to destroy the bacteria. This direct correlation between the amount of bacteria in the wound and the neutrophil concentration also affects the reduction of neutrophils from the chronic wound. The fraction (the last item on the right side of (2.5)), which is also mentioned in ^[26], describes this process. Moreover, two bounded smooth functions $G_N(O(x, t))$ and $H_N(O(x, t))$ are defined as

and

In according with (2.3), the boundary conditions for (2.5) are

The neutrophils are assumed to be originally concentrated heavily near the edge of the chronic wound, the initial condition for (2.5) is

where $\epsilon$ is the distribution coefficient.

2.3.   Dimensionless bacterial concentration model

The movement process of the dimensionless bacterial concentration in a chronic wound is governed by the following PDE model ^[26]

where $D_B$ is the diffusion coefficient of bacteria; $\alpha_{B}$ is the proliferation coefficient of bacteria; $K_O$ is the influence coefficient of oxygen; $\theta$ is the influence coefficient of other leukocytes; $K_{NR}$ is the influence coefficient of neutrophils; $\lambda_{RB}$ is the attenuation coefficient of neutrophils' influence; $\lambda_B$ is the natural death rate of bacteria. Bacteria are removed by neutrophils. An increase in the neutrophil concentration within the wound site will lead to the destruction of bacteria. Thus, $K_{NR}$ denotes the influence coefficient of neutrophils for bacterial destruction while $\lambda_{RB}$ denotes the attenuation coefficient of neutrophils' influence caused by bacteria. Since the bacterial concentration satisfies the natural boundary condition (a zero flux condition) at $x = 0$ and $x = 1$, the boundary conditions for (2.10) are

Meanwhile, the bacteria are assumed to be initially concentrated densely in the center of the chronic wound within $[0, \epsilon]$ and nonexistent elsewhere. It is corresponding to the neutrophils, which are assumed to be initially concentrated heavily near the edge of the chronic wound. Then, neutrophils are chemically attracted by the chemoattractant to destroy bacteria. Moreover, the assumption that bacteria are initially concentrated in the center of the chronic wound is mentioned in ^[26]. In practice, the process of bacteria slowly moving from a central area of a wound to the outer edge is reasonable. Thus, the initial condition for (2.10) is

where $\epsilon$ is the distribution coefficient like in (2.9).

2.4.   Dimensionless chemoattractant concentration model

The movement process of the dimensionless chemoattractant concentration in a chronic wound is governed by ^[27]

where $D_C$ is the diffusion coefficient of chemoattractant; $\alpha_{C}$ is the proliferation coefficient of chemoattractant produced by the bacteria; $\lambda_C$ is natural death rate of chemoattractant. In accordance with (2.11), the boundary conditions for (2.13) are

In accordance with (2.12), the initial condition for (2.13) is

where $\epsilon$ is the distribution coefficient like in (2.9).

2.5.   The optimal control problem

Since bacterial concentration causes a serious threat to health cure, the hyperbaric oxygen therapy $u(t)$ must be manipulated carefully to inhibit the growth of bacterial concentration and even eliminate bacteria at fixed treatment time $t = T$. Gill et al. mention that hyperbaric oxygen therapy is expensive. The cost of this therapy depends on the amount of time it takes to use oxygen. We hope to achieve better therapeutic results at a lower economic cost, which has important socioeconomic significance. The cost functional is designed to minimize the bacterial infection as well as the amount of time spent in the hyperbaric oxygen therapy. To this end, we consider the following cost functional

where $\omega$ is the weight coefficient, which depends on the treatment cost. The first term of cost functional (2.16) denotes the bacterial concentration at fixed treatment time $t = T$ and the second term denotes the economic cost by hyperbaric oxygen therapy. This paper considers the inhibition of bacterial propagation and even completes the elimination of bacteria to achieve a smooth healing of a chronic wound through the hyperbaric oxygen therapy. Thus, bacterial concentration is the most direct variable and an objective evaluation criterion for therapeutic effects. Other variables, including oxygen, neutrophils and chemoattractants, are coupled with bacteria through PDEs. It is difficult to obtain the direct quantitative relationship between bacteria and indirect variables and then hard to include these indirect variables in the cost functional to obtain an objective evaluation criterion. For example, chemoattractant attracts neutrophils to destroy bacteria while chemoattractant is produced by bacteria. Thus, the chemoattractant concentration is difficult to be quantified in the cost functional. Moreover, the proposed objective function is also used in ^[21].

Thus, we obtain the following PDE optimization problem: given the dimensionless oxygen concentration model (2.1), dimensionless neutrophil concentration model (2.5), dimensionless bacterial concentration model (2.10) and dimensionless chemoattractant concentration model (2.13) with corresponding boundary conditions and initial conditions, choose the control $u(t)$, in accordance with (2.2), to minimize the cost functional (2.16). We refer to this problem as Problem P$_0$.

3.   The method of lines

Problem P$_0$ is a PDE optimization problem. In order to simplify the optimal control design, we apply the method of lines, which only discretizes the space while keeps the time continuous, to reduce the PDEs (2.1), (2.5), (2.10) and (2.13) into ODEs ^[28]. First, we decompose the whole space $[0, 1]$ into equally-spaced intervals $x_i-x_{i-1} = 1/M$, $i = 1, ..., M$, where $M$ is an integer of spatial discretization, $x_0 = 0, \; x_M = 1$ and define

Then, based on the (3.1), we obtain the following finite difference approximations by center difference

and

Thus, based on the (3.2) and (3.3), the dimensionless oxygen concentration model (2.1), dimensionless neutrophil concentration model (2.5), dimensionless bacterial concentration model (2.10) and dimensionless chemoattractant concentration model (2.13) become

Under the spatial discretization (3.1), the boundary conditions (2.3) and (2.8) become

The boundary conditions (2.11) and (2.14) become

The initial conditions (2.4) and (2.9) become

The initial conditions (2.12) and (2.15) become

Moreover, the cost functional (2.16) becomes

Now, we state our approximate problem as follows: given the ODEs (3.4) with the initial conditions (3.7) and (3.8), choose the control $u(t)$, in accordance with (2.2), to minimize the cost functional (3.9). We refer to this problem as Problem P$_1$.

4.   Optimal switching time control

4.1.   Time-scaling transformation approach

Problem P$_1$ is an ODE optimization problem, where the control $u(t)$ is limited to 0 or 1. In the hyperbaric oxygen therapy, we should choose the optimal time intervals applied by the oxygen. Note that the time of oxygen inhalation should not be too long, no more than three times a day, no more than 120 minutes (2 hours) each time, otherwise, it is easy to cause oxygen toxicity for a human ^[29]. First, the whole time horizon $[0, T]$ is divided into $\gamma$ time intervals ^[30]

where $t_k, \; k = 1, \dots, \gamma-1,$ are time switching points. However, unequal time intervals $[t_{k-1}, t_{k}]$, $k = 1, ..., \gamma$ make the ODE optimization problem difficult to deal with, where the optimized decision variables are the time intervals. Thus, we apply the time-scaling transformation approach to transform the ODE optimization problem with changed switching time into a new optimization problem with fixed switching time ^[31,32]. The relationship between $t$ (changed switching time) and $s$ (fixed switching time) shown in Figure 2 can be also defined as

where $t(0) = 0$ and $t(\gamma) = T$; $\zeta ^k$, $k = 1, ..., \gamma$, are time interval variables, which should be optimized. Let $\boldsymbol \zeta = [\zeta ^1, \zeta ^2, ..., \zeta ^\gamma]^\top$, where $\zeta ^k = t_k - t_{k - 1}, k = 1, ..., \gamma$ and

and when $u(t) = 1$, $0 \leq\zeta ^k\leq T_h,$ where $T_h$ denotes 2 hours.

Integrating (4.2) over $[0, s]$, we can obtain the time-scaling function as

where $\left\lfloor s \right\rfloor$ donates an integer which is not larger than $s$. Accordingly, we should change the Problem P$_1$ with time variable "$t$" into a new optimization problem with auxiliary variable "$s$". Let denote

Combined with the chain rule, the reduced dimensionless oxygen concentration model (3.4) becomes

Based on the (4.6), (3.4) becomes

where $s \in (k-1, k), k = 1, \dots, \gamma.$ The initial conditions (3.7) and (3.8) become

After the time-scaling transformation approach (4.2), the cost functional (3.9) becomes

Finally, the optimal parameter selection problem is stated as: given the ODEs (4.7) with the initial conditions (4.8), choose the optimal switching time points $t_k, \; k = 1, \dots, \gamma-1$ combined with the control $u(t)$ to minimize the cost functional (4.9). We refer to this problem as Problem P$_\gamma$.

4.2.   Solving problem P$_\gamma$

Problem P$_\gamma$ is a nonlinear programming problem. Thus, we should compute the gradient of the cost functional (4.9) corresponding to the time interval variables $\zeta^k, k = 1, ..., \gamma$. However, computing its gradient is a non-trivial task since (4.9) is an implicit, rather than an explicit function.

For each $\mu = 1, 2, \ldots, \gamma$ and $i = 1, ..., M-1$, we apply the integration and obtain

If $k\leq \mu$, we differentiate (4.10) corresponding to $\zeta_k$ and obtain

where $\delta_{k\mu}$ denotes the Kronecker delta function. If $k > \mu$, $\zeta_k$ does not affect the $\tilde O_i(s)$ during $s\in[{\mu-1}, \mu)$, thus we have

Then, we differentiate (4.11) corresponding to $s$ and obtain

From (4.10)–(4.13), we can obtain the state variations $\frac{ {\partial \tilde O_i(s)}}{\partial \zeta_k}$ for each $k = 1, 2, \ldots, \gamma$, $i = 1, ..., M-1$

where $s\in[{\mu-1}, \mu),$ $\mu = k, k+1, ..., \gamma$.

For the initial condition (4.8), we have

Similarly, we can obtain the state variations $\frac{ {\partial \tilde N_i(s)}}{\partial \zeta_k}$, $\frac{ {\partial \tilde B_i(s)}}{\partial \zeta_k}$ and $\frac{ {\partial \tilde C_i(s)}}{\partial \zeta_k}, k = 1, 2, \ldots, \gamma$, $i = 1, ..., M-1$.

Moreover, for the cost functional (4.9), the gradient formulas corresponding to $\zeta_k, k = 1, 2, \ldots, \gamma$ can be obtained using the chain rule of differentiation

Based on the gradient formulas (4.16), the existing nonlinear programming algorithm, such as Sequential Quadratic Program (SQP) algorithm ^[33], can be applied to solve the Problem P$_\gamma$. The flow chart is shown in Figure 3.

5.   Computational numerical analysis

For the computational numerical analysis, the proposed method in Section 4 is implemented in MATLAB software, which uses ODE15s to solve the reduced ODEs (4.7) for the cost functional and gradient formulas and then apply the SQP algorithm to perform the optimization steps. The parameter values for the dimensionless oxygen concentration model and dimensionless neutrophil concentration model are set as: $\eta = 0.2284$, $\kappa = 5.4$, $\lambda_{NO} = 37$, $\lambda_{BO} = 22.7872$, $\lambda_{O} = 2.4667$, $D_N = 0.02$, $\kappa_N = 10$, $\alpha_{BN} = 14.28$, $\lambda_N = 5$ and $e = 30$. Meanwhile, the parameter values for dimensionless bacterial concentration model and dimensionless chemoattractant concentration model are set as: $D_B = 0.0001$, $\alpha_{B} = 1.26$, $K_O = 0.75$, $K_{NR} = 2$, $\lambda_{RB} = 3.73$, $\lambda_B = 5$, $D_C = 1.5$, $\alpha_{C} = 20$ and $\lambda_C = 0.9$, which are consistent with ^[21]. Moreover, the distribution coefficient $\epsilon$ is set as $0.01$ and the whole time horizon $T$ is set as $0.432$, which corresponds to $24$ hours.

The spatial discretization number $N$ is set as $61$ and the weight coefficient $\omega$ is set as $0.07$. Our computational numerical analysis is carried out on a personal computer with the following configuration: Intel(R) Core(TM) i7-6820HQ CPU, 2.70GHz CPU, 16.00GB RAM, 64-bit Windows 10 Operating System.

We divide the whole time horizon into $\gamma = 5$ time intervals and set corresponding control sequences $u(s)$ as $0$, $1$, $0$, $1$, $0$. The initial guess for $\zeta ^k, k = 1, ..., \gamma$ is $T/\gamma$, which denotes the equal time intervals. The optimal switching time control sequences under $\gamma = 5$ are shown in Figure 4a and the dimensionless bacterial concentration at $x = 0$ under $\gamma = 5$ is shown in Figure 4b. Hyperbaric oxygen therapy is not required due to the high concentration of oxygen at the initial time. Then, it is applied when the oxygen concentration is reduced but the bacterial concentration is still relatively high. Finally, our normal human immune system completely replaces expensive treatments to inhibit or even destroy these bacteria when the bacterial concentration is low. The distribution of the dimensionless oxygen concentration, dimensionless neutrophil concentration, dimensionless bacterial concentration and dimensionless chemoattractant concentration is shown in Figure 5a, 5b, 5c and 5d, respectively. The oxygen concentration (shown in Figure 5a) decreases first, then increases and then decreases due to the optimal switching time control sequences (shown in Figure 4a). The neutrophils (shown in Figure 5b) remain dense near the edge of the chronic wound and are chemically attracted to the center of the chronic wound very slowly. Since the chemoattractant (shown in Figure 5d) does not reach the edge of the chronic wound, neutrophils are lack of movement. Through the regulation of oxygen, the bacterial concentration approaches to 0. Correspondingly, the chemoattractant concentration drops to 0 due to the quick loss of the bacteria.

We extend the number of time intervals into $\gamma = 7$ and set corresponding control sequences $u(s)$ as $0$, $1$, $0$, $1$, $0$, $1$, $0$. The optimal switching time control sequences are shown in Figure 6a and the dimensionless bacterial concentration at $x = 0$ is shown in Figure 6b. The distribution of the dimensionless oxygen concentration and dimensionless bacterial concentration is shown in Figure 7a and 7b, respectively. Since the optimal switching time control sequences ($\gamma = 7$) are similar to the case of $\gamma = 5$, the dynamics of oxygen, neutrophils, bacteria and chemoattractants are almost identical to those of $\gamma = 5$. However, note that increasing $\gamma$ from $5$ to $7$ does not lead to any significant change in the inhibition of bacteria. Thus, choosing $\gamma = 5$ for the time intervals is sufficient.

The results also provide protocols for switching treatment on/off in practice and their meaning in terms of bacterial suppression. From Figures 5a, 5c, 7a and 7b, we can see that when the dimensionless oxygen concentration is less than about 0.75 and the dimensionless bacterial concentration is more than about 0.45, the hyperbaric oxygen therapy is applied effectively. Otherwise, our normal human immune system can inhibit or even destroy these bacteria and the hyperbaric oxygen therapy is not necessary. Note that there are much uncertainty and complex coupling relationship between variables in practice, which will influence the optimal therapeutic effect. Moreover, ^[31] presents the convergence proof that the optimized control sequences after the time-scaling transformation can approach the theoretic optimal control as $\gamma\rightarrow +\infty$. When the number of time intervals is increased to a certain extent, the optimal results are similar and not sensitive to $\gamma$. However, we cannot choose $\gamma$ too large due to consideration of human safety and the limitations of medical conditions. We choose $\gamma = 5$ (twice a day) and $\gamma = 7$ (three times a day) in the computational numerical analysis. Our numerical results show that the optimal switching time control sequences are similar when $\gamma$ is increased from $5$ to $7$, which means the results are not very sensitive to $\gamma$.

Different from traditional empirical strategies, this paper's work can accurately give better protocols for switching treatment on/off at a lower economic cost. Through precise hyperbaric oxygen regulation, bacterial propagation is effectively inhibited, thereby the wound healing is improved. The results of the proposed method can provide some guidance for the medical staff by use of the hyperbaric oxygen therapy for wound healing in practice. Moreover, combined with image recognition technology, our proposed method can be easily extended to the intelligent medical field.

6.   Conclusion and future work

In this paper, we propose the optimal switching time control of the hyperbaric oxygen therapy for bacterial concentration suppression. The bacterial concentration is coupled with the oxygen concentration, neutrophil concentration and chemoattractant concentration governed by PDEs. The method of lines is first used to reduce the PDE optimization problem into an ODE optimization problem. Then, the time-scaling transformation approach is applied and the gradient formulas of the cost functional corresponding to the time intervals are derived based on the sensitivity analysis. Finally, a new effective computational method is proposed and computational numerical analysis shows the effectiveness of the proposed strategy. In the future, we will work on the experiments to further verify the feasibility of the proposed strategy. Moreover, the iterative learning control ^[34], control parameterization with optimized control values ^[35,36] and adaptive control ^[37] for the hyperbaric oxygen therapy can be considered.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China 61473253, the Medical Science and Technology Project, Zhejiang Province (2016154240, 2017192644).

Conflict of interest

All authors declare no conflicts of interest in this paper.