Effective Spatio-Temporal Regimes for Wound Treatment by Way of Macrophage Polarization: A Mathematical Model

Wound healing consists of a sequence of biological processes often grouped into different stages. Interventions applied to accelerate normal wound healing must take into consideration timing with respect to wound healing stages in order to maximize treatment effectiveness. Macrophage polarization from M1 to M2 represents a transition from the inflammatory to the proliferative stage of wound healing. Accelerating this transition may be an effective way to accelerate wound healing; however, it must be induced at the appropriate time. We search for an optimal spatio-temporal regime to apply wound healing treatment in a mathematical model of wound healing. In this work we show that to maximize effectiveness, treatment must not be applied too early or too late with respect to peak inflammation. We also show that the effective spatial distribution of treatment depends on the heterogeneity of the wound surface. In conclusion, this research provides a possible optimal regime of therapy that focuses on macrophage activity and a hypothesis of treatment outcome to be tested in future experiments. Finding optimal regimes for treatment application is a first step toward the development of intelligent algorithms for wound treatment that minimize healing time.

Wound healing consists of a sequence of biological processes often grouped into different stages. Interventions applied to accelerate normal wound healing must take into consideration timing with respect to wound healing stages in order to maximize treatment effectiveness. Macrophage polarization from M1 to M2 represents a transition from the inflammatory to the proliferative stage of wound healing. Accelerating this transition may be an effective way to accelerate wound healing; however, it must be induced at the appropriate time. We search for an optimal spatio-temporal regime to apply wound healing treatment in a mathematical model of wound healing. In this work we show that to maximize effectiveness, treatment must not be applied too early or too late with respect to peak inflammation. We also show that the effective spatial distribution of treatment depends on the heterogeneity of the wound surface. In conclusion, this research provides a possible optimal regime of therapy that focuses on macrophage activity and a hypothesis of treatment outcome to be tested in future experiments. Finding optimal regimes for treatment application is a first step toward the development of intelligent algorithms for wound treatment that minimize healing time.

INTRODUCTION
Delayed wound healing presents an important health-care problem. There is no decisive finding regarding the best therapy for delayed wound healing due to the variety of complications that can ensue [1]. In the case of acute wounds, a growing research area in wound healing has focused on methods to accelerate wound healing such as application of an electric field, application of stem cells, and the passive release of therapeutic molecules in so-called smart bandages [2]. Less attention has been given to the timing of any given therapy. A treatment may exert no effect or even a negative effect on healing tissues if not applied appropriately [3]. For example, certain treatments can accelerate specific stages of inflammation but have little effect on others [4,5]. In other cases, treatment can induce toxic side-effects [6]. For this reason each treatment should be applied only at the appropriate stages of wound healing. The emergence of bioelectronic devices provides the opportunity to achieve drug delivery in a continuous and controlled fashion [7,8].
Wound healing consists of several stages: hemostasis, inflammation, proliferation, and remodeling [9][10][11]. During hemostasis, a blood clot is formed to stop bleeding [12]. Inflammation begins immediately after a blood clot formation. At this time, chemokines are released to attract immune cells -neutrophils and macrophages. During the proliferative stage, angiogenesis, collagen deposition, and the formation of granulation tissue occur. The main participants in this stage are fibroblasts and anti-inflammatory macrophages. During remodeling, unnecessary blood vessels are removed, the extracellular matrix is remodeled, and tissue architecture is restored.
Wound acceleration can be achieved by shortening the duration of one or more of the wound healing stages. Determining the optimal timing for a therapy requires one to map the direct effect of the therapy on the targeted biological processes to the overall wound closure time.
Mathematical modeling in this field may help to predict optimal treatment regimes and to plan future experiments. Existing mathematical models have served to investigate different aspects/stages of wound healing [13][14][15]. Models of inflammation have suggested strategies to avoid chronic inflammation by control of neutrophil apoptosis and macrophage phagocytosis [16]. The model presented in Xue et al. [17] suggests a mechanism by which a deficiency in oxygen supply can limit macrophage recruitment and slow healing. The models describing production of the extracellular matrix by fibroblasts allowed researchers to investigate scar formation [18]. The role of cell migration and proliferation on wound closure was investigated in Javierre et al. [19]. Thus, we propose that a qualitative model with a proper level of abstraction can be used to predict the response of a wound to dynamic therapy.
In this study, we investigate mathematically the most effective spatio-temporal regimes for drug delivery to accelerate wound closure. In particular, we identify the effects of delivering substances that can accelerate macrophage polarization on wound closure times. Macrophage polarization modeling, to our knowledge, has been applied to several biological situations but not to wound healing [20,21].
Macrophages play important roles at all stages of wound healing [9][10][11]22]-from clearing the wound of debris by phagocytosis to maintaining cell proliferation of the tissue being repaired. Macrophages can perform various tasks due to their ability to transform into several phenotypes depending on external stimuli [23]. The most famous phenotypes are M1 (proinflammatory) and M2 (anti-inflammatory), identified in vitro [24]. In a wound, the set of stimuli received by macrophages is constantly changing and their phenotype undergoes a dynamic transition. The ability to regulate the macrophage phenotype and achieve a controlled time-dependent transition from the M1 phenotype to the M2 phenotype is a promising approach to accelerate wound healing [25][26][27]. There are several substances that can induce macrophage polarization from the pro-inflammatory to the anti-inflammatory subtype [28][29][30]. This type of treatment can have adverse effects if applied improperly [31,32]. For example, a fast transition from the inflammatory to the proliferative stage may slow down the cleaning of the wound of debris (e.g., removal of harmful bacteria and damaged cells).
We focus our attention on accelerating wound closure of acute wounds. For simplicity, and to this end, we consider a simple model with only one stable state -the healthy one. That is, our model does not capture switching between chronic and normal wound healing regimes. However, in order to appropriately capture trade-offs of early treatment, we model wound debris over time, which is actively degraded by M1 macrophages. Tracking this state, we quantify wound cleaning time with the understanding that any remaining wound debris can be an indicator of prolonged inflammation and potential infection preventing wound closure. That is, wound debris should decay in a timely manner for the treatment to be realizable.
In summary, we examine wound healing time and wound debris cleaning time in response to different spatio-temporal signals inducing macrophage polarization. Overall, decoupling the different modalities of wound healing trajectories reduces complexity of the model and allows us to gain intuition for optimal treatment strategies. We find that actuation of M1-M2 polarization must be applied with care and optimal timing can depend on the duration of the treatment, time of initiation, placement of the actuator, and initial distribution of wound debris.

Background
The timing and coordination of biological processes involved in wound healing are complex, however, many of them can be lumped into a single state [33]. A good model should contain the minimal number of variables needed to describe the process [14].
We consider a five state system, which provides an appropriate level of abstraction for our study. In order to study the effects of macrophage polarization we must make them explicit states. M1 macrophages are associated with the early stages of inflammation [11,34] and the upregulation of M2 signals corresponds to initiation of the proliferation stage [9,10]. The remaining states defined as wound debris, temporary tissue, and new tissue represent primary activity in early and later stages of wound healing.
Debris of both internal and external origin appear in the wound at the moment of injury, when healthy tissue is damaged. The debris activate M1 macrophages and are then removed by them.
The proliferative stage is characterized by the production of an extracellular matrix facilitated by fibroblasts and the growth of new blood vessels [39,40]. These processes are regulated by M2 macrophages [10] and are needed to help new tissue form FIGURE 1 | (A) Geometry of the model: the r-axis is directed from the wound center; R is the wound radius. (B) Schematic of the model of wound healing. Wound debris D b attract M1 macrophages that remove debris. M1 macrophages become M2 macrophages, which induce production of temporary tissue C that helps new tissue N to grow.
correctly, thus, we consider these agents as "temporary tissue" in our model.
During the remodeling stage, the extracellular matrix is remodeled and extra blood vessels are removed [39]; new healthy tissue replaces temporary tissue indicating wound closure. Here, we consider a wound where only the top epithelial layer is damaged and new tissue grows from the edge of the wound as epithelial cells proliferate and migrate inward. Thus, new tissue in our model is growing as a sheet from the wound edge [41].
These five states capture the well characterized wound healing stages and changing wound conditions inducing a biological response.

Equations
Consider a wound of radius R. Let the axis r to be directed from the wound center to the edge ( Figure 1A). Concentrations of substances and populations of cells are functions of r. A schematic of the interaction between the modeled biological processes including macrophage participation in wound healing is shown in Figure 1B. D b is wound debris consisting of damaged cells and bacterial cells promoting infection. We assume the wound debris to be non-active and only eliminated by M1 macrophages: where M 1 is the population of M1 macrophages, which are attracted by debris (k 2 D b ) and removed in the reactions corresponding to debris elimination ( ), and natural death (k d1 M 1 ). Spatial migration of macrophages is described by a classical diffusion term [42]. Thus, the dynamics of M1 macrophages are described by the following equation: where the diffusion is written in a cylindrical coordinate system. M 2 is the population of M2 macrophages whose dynamics are driven by M1 polarization (k 4 and migration as follows: Proliferation is a complex process involving fibroblasts. The extracellular matrix (ECM) formation and its partial destruction is controlled by a complex coordination of enzymes [39,40]. This finely tuned system supporting new tissue growth has been modeled in several works [18,43,44]. Here, we use a simplified model of proliferation. The temporary tissue variable C represents the temporarily formed ECM and active enzymes. Temporary tissue production is induced by M2 macrophages (k 5 M 2 ) and is destroyed thereafter (−k r C). The dynamics are described by the following equation: The state C is an intermediate state leading to the growth of new healthy tissue, N. To model growth of new healthy tissue, we assume that the top layer of the skin, the epithelium, grows as a sheet from the edge of the wound [41]. Epithelial cells divide, but this process slows down if there are too many cells. Thus, the epithelial growth velocity is described as N(1 − N). If the amount of new tissue reaches N = 1 it does not grow any more. New tissue cells can migrate to adjacent areas with a migration rate D n that is much slower than that of macrophages. The equatioṅ N = aN(1 − N) + ∂ 2 N ∂x 2 is known in the theory of nonlinear dynamical systems as the Fisher-Kolmogorov equation [45], and its solution can be represented as a running wave [42]. However, this is possible only in the presence of temporary tissue, so the rate of new tissue growth is proportional to C. The dynamics of new tissue are described by the following equation: In order to reduce the number of parameters, we reparameterize the model by introducing new variables: where T and L are characteristic time and length scales. The system of equations in the reparameterized form is: where The parameter values used in numerical simulations are listed in Table 1.
We assume a time scale of 8 h and a spatial scale of 0.03 mm. For β the rate of attraction of macrophages by debris/pathogen was reported in Hau et al. [46] as 1/day and in Andersson et al. [47] as 0.5/day. We set β = 1. The degradation rate of macrophages varies in literature from 0.014/h [34] to 0.12/h [49]. We accept the values γ 1 = γ 2 = 0.1. We select parameters k, q and ρ such that M1 and M2 dynamics agree with those found in Du et al. [48]. The value of µ was set such that temporary tissue is removed 3 weeks after injury [39]. The coefficient of diffusion for macrophages can be deduced from cell migration experiments. In Wheeler et al. [51], the mean displacement of cells on plastic was 20-30 µm in 6.5 h, giving the estimate for the diffusion coefficient [50]D=0.27 − 0.61. We assumeD = 0.32.
Finally, the parametersD n = 0.0003 andα = 1.8 were selected in such a way that complete wound closure was observed by day 14, while no changes in wound closure were observed in the first 2 days as is consistent with some experiments [52]. The sensitivity analysis of the model to variations in parameters is shown in Supplementary Figures S5, S6.
The initial conditions are: We assume zero-flux boundary conditions for macrophages on the right and left boundaries of the considered region: New tissue is assumed to be constant at the edge of the wound and non-moving through the center of the wound: ∂n ∂r |r =0 = 0, n|r =R/L = 1.
The system of Equations (6)-(10) was solved numerically on a uniform mesh consisting of 100 spatial cells. Five equations were written in each cell, with the diffusion term approximated using a central difference scheme. The resulting system of 500 ordinary differential equations was solved in Matlab R2020a by the ode15s solver. The results of wound healing model simulations are shown in Figure 2. The wound radius is measured as the distance from wound center to the location r, where n(r) > 0.95.

Model of Wound With Actuator
In order to investigate regimes of wound healing treatment we include actuator induced macrophage polarization into the model. This actuator is applied at a radius r = r p . We assume that the actuator delivers a biochemical at the point of application and its concentration θ is described by the following function in time and space (see Figure 3A): where θ 0 is the amplitude of the treatment controlled by the actuator. One can see that θ = θ 0 at the location of the actuator,  r = r p . In this work we assumed σ 1 = 0.3 mm and σ 2 = 0.09 mm. The treatment substance delivered affects macrophage polarization, so equations for m 1 and m 2 with actuators may be rewritten as follows: In order to find the optimal regime, we test the model with impulses of actuator treatment of duration t beginning at time t 0 (see Figure 3B): We present a piece-wise linear function for simplicity. We assume linear growth and a decrease at the beginning and end of the signal. We tested other impulse-like shapes for actuation and found that the general behavior of the system response remained unchanged.

RESULTS
We define wound healing time as the time from injury (t=0) to the moment when the wound radius reaches zero. Application of an actuator that accelerates macrophage polarization in the model decreases the time of wound healing. The results are shown in Figures 4-6.
The beginning time of actuation plays an important role in wound healing (see Figure 4). Large values of t 0 make the treatment less effective: wound healing time increases as t 0 increases. In our simulations for t 0 > 3d the treatment does not have any effect: the value of healing time tends to the value of healing without treatment (13.47 days for the given set of parameters).  Interestingly, for shorter treatment durations t, healing time plots have a minimum away from the extremes (see plots for t < 2d in Figure 4A). There exists a timewindow in the wound healing process when treatment is most effective. This time-window is likely characterized by some underlying biological process. If the treatment duration is short and applied early, the treatment ends before the most effective time-window is reached. This is not observed for longer treatments ( t = 2 − 3 days) because even with early application at t 0 = 0, the treatment duration overlaps with the most effective time-window for treatment.
This means that there is a non-trivial optimal treatment beginning time t 0 . In other words, there is a short window of time during wound healing when artificial acceleration of macrophage polarization is most effective. The plots for other actuator positions r p are shown in the Supplementary Figure S1. Similar trends are observed for different placement of the actuators.  Figure S2. Again, we find a similar trend regardless of actuator placement. Figure 6A shows the dependence of wound healing time on the actuator position r p . The treatment substance in this model is approximated by a piece-wise function (14) on a bounded domain. To ensure the bounded domain stays within the region of interest we constrain the actuator placement: 0.3 mm≤ r ≤ 2.7 mm. The plots are decreasing from the wound center to the edge. This implies that an actuator located close to the wound edge is beneficial. The plots for other t and t 0 are shown in Supplementary Figures S3, S4.
The results above might be the consequence of the initial uniform distribution of the debris in the wound (see initial conditions). We provide additional simulations with debris accumulation in the center and at the edge of the wound. Alternative initial conditions for the debris variable take the following form: and a(r)| t=0 = 2r/R.
Wound healing time dependence on r p for the three different types of initial conditions on debris distribution is shown in Figure 6B. One can see that debris distribution in the wound bed affects the dependence of healing time on actuator placement. For the case a(r) = 2r/R, there is more debris on the edge. If the actuator is placed near the edge, we can get much shorter healing times. For the case 2(1 − r/R) (more debris in the wound center), there exists an optimal position of the actuator away from the edge (r p ≈ 1.8 mm). One can see that the optimal actuator position is sensitive to the distribution of the debris in the wound. Debris distribution cannot be estimated in the framework of this rough modeling and must be investigated experimentally. Figure 7 demonstrates the limitations of treatment regimes. In addition to wound healing time, we define the wound cleaning time for debris removal. Because debris is a diminishing variable in our model, it tends to zero as time goes to infinity. We define wound cleaning time as the time when the maximal value of debris across the wound bed falls below a small threshold 0.025: max r a < 0.025. (20) Because the treatment applied in this model accelerates M1 to M2 transition, it not only accelerates the proliferation-remodeling stages but removes M1 cells performing debris removal. Too large application of this treatment can make debris removal too slow-this is the cost of accelerating the proliferative and remodeling stages. Figure 7A shows that a shorter wound healing time corresponds to a longer time required for wound cleaning or debris removal. For some scenarios, wound cleaning may take even longer than wound closure, which is physically unrealizable and indicative of complications in wound healing.
We can see two main limitations of the treatment regimes, divided by vertical dotted lines in Figure 7A. The regimes with small t 0 (left of the first dotted line) lead to too slow wound cleaning. Regimes with very large t 0 (to the right of the second dotted line) result in insufficient improvement in wound healing. We define "insufficient" by a less that 10% reduction in wound healing time when treatment is applied. This threshold is of course arbitrary and can be chosen by the user. For reference, in Liang et al. [53], the authors demonstrated accelerated wound closure on a pig wound by applying a continuous external electric field and reduced the time for wound closure by ∼12-18%. Thus, we opine that below a 10% threshold, the cost and potential side effects of applying the treatment outweigh the benefit.
Only the regimes with middle values of t 0 , between the vertical dotted lines in Figure 7A, shorten wound healing time sufficiently, while keeping a reasonable wound cleaning time.

Figure 7B
shows these 3 types of regimes in the (t 0 , t) parametric plane. The regimes marked as red squares correspond to regimes when wound cleaning takes a longer time than wound healing, whereas blue diamonds correspond to regimes when wound healing time diminishes less than 10% in comparison with a non-treated wound. Green crosses represent the effective regimes of wound treatment. We find that the treatment should be applied between 0 and 1 days. Applying treatment past 1 day in all scenarios does not significantly reduce the wound healing time, hence, the treatment is ineffective. We note that in our simulations, M1 macrophage activity, which is associated with inflammation, peaks at around 24 h. Thus, we may hypothesize that treatment is best applied shortly after or just shy of peak inflammation, thereby, accelerating the transition from the inflammatory state to the proliferative state. Applying treatment too early can result in chronic wounds given a less than reasonable time is provided for inflammation, a stage in wound healing critical for preventing infections.
Of course, one can choose more stringent conditions for optimal regimes (maybe 10% healing time improvement is not enough, or the threshold in the condition (20) should be smaller. Some variants of other threshold selections are shown in Supplementary Figure S7). Among all the regimes presented in Figure 7B the smallest healing time (11.37 days) is observed for the treatment regime t 0 = 8 h and t = 2 days 8 h. This reduces wound healing time by approximately 2 days. In practical situations there might be additional constraints on the duration of treatment, or maximum concentration that is not considered here. However, our model demonstrates the principles of wound treatment regime optimization.

DISCUSSION
To our knowledge, the mathematical model presented here is the first one that considers the role of macrophage polarization explicitly in wound healing. The model takes into account the presence of two types of macrophages in the wound-M1 and M2. It is believed that in addition to these primary red squares -wound cleaning takes a longer time than wound healing; blue diamonds -wound healing time diminishes less than 10% in comparison with the non-treated wound (r p =0.6 mm). Arrows labeled "a" and "b" represent the same change in treatment strategy but demonstrate different outcomes due to different reference points.
subtypes, there are several other subtypes of macrophages. All subtypes of macrophages were found in vitro as a result of identifying corresponding activating stimuli. The exact subtypes of macrophages in wounds have not been established and are considered to be roughly similar to those found in vitro. It is most likely that macrophages of different types can be present in the wound simultaneously. However, the functions that these macrophages perform are more important for the healing process than the markers found in vitro. Therefore, in this model, we clearly divide the functions of the macrophages into inflammatory and reparative, keeping in mind that such pure cell lines may not exist in reality. This separation of functions of macrophages gave us the opportunity to draw up a rather simple naive model and identify general patterns for the effect of varying treatment regimens on wound healing time.
Admittedly, this model has its limitations. The model constructed in this work does not take into account the details of macrophages polarization due to the poor knowledge of this mechanism in vivo [34]. It is known that M2 macrophages can appear in response to stimulation by certain cytokines, for example, IL4, produced by basophils and mast cells [11]. There are indications that M1 macrophages are converted to M2 subtypes after they phagocyte apoptotic neutrophils [38]. We made a model in which the M1 population is replaced by M2 although the mechanisms of this transition in vivo are not well understood. Thus the mechanisms of polarization and their effects on timing cannot be included. Still we derive biologically meaningful results.
Our results imply that a treatment targeting macrophage polarization should not take place too early. Otherwise, M1 macrophages do not have enough time to eliminate wound debris and wound healing can be delayed. Details of this scenario is out of scope of our model. The onset of this scenario is especially important for infected wounds when the debris consists not only of damaged cells but of bacterial pathogen [54]. In the case of bacterial pathogen, the first equation in the model should include pathogen reproduction and regimes where pathogen persistence may occur. This may lead to continuation of M1 recruitment, inflammation persistence, and prevent full wound closure. On the other hand, treatment targeting macrophage polarization should not begin too late, because it becomes ineffective and toxic side effects are unknown.
The model can still be used to guide experiments. The parametric plane of treatment regimes in the (t 0 , t) coordinate plane ( Figure 7B) may serve as a first approximation for planning wound treatment strategies. Although the boundaries separating the three regimes in Figure 7B may shift across experiments, we might expect the general pattern to hold. For example, if an experimentally tested treatment regimen is found to not shorten healing time to a desirable extent, then the implication is that the domain of actuation is in the blue region of Figure 7B. A reasonable choice for the next strategy is to decrease t 0 and increase t with the goal of shifting the system from the blue to green region (arrow a in Figure 7B). However, if the new regimen leads to complications in the inflammatory stage, this implies the shift pushed the actuation domain into the red region, as shown by arrow b in Figure 7B. Having this mathematically derived parametric plane may help experimenters orient in an unknown field of regimes. The specific values of the thresholds and the corresponding shapes of the boundaries can be clarified experimentally.
We note that the geometry of the domain of actuation corresponds to a ring of a given width at a given radius. Recent advances in the design of bioelectronic devices provide the ability to deliver biochemicals with spatio-temporal precision [7]. Thus, we would like to take full advantage of these capabilities to understand the dependence of treatment effectiveness on proximity to the wound center compared to the wound edge. The lack of knowledge about debris distribution in the wound makes it difficult to predict the optimal actuator position. However, we demonstrate several possible solutions for different debris distributions in the wound. Experiments on wound treatment should be done to clarify the optimal spatial distribution of treatment actuators.
Nonetheless, the general trends for changes in wound healing times with respect to changes in actuation timing are similar for various fixed spatial profiles. The importance of considering spatial context in wound treatment remains an open question. If considering spatial information can significantly improve treatment strategies, then such numerical studies can inform future actuator design.

CONCLUSION
Wound healing consists of a sequence of stages with different cells performing different functions. Therefore, it is reasonable to assume that at each stage, different medications should be applied to improve healing. However, to our knowledge, there are not many studies of dynamic wound treatment regimens. In this work, we attempted to find the optimal wound treatment regimen by affecting macrophage polarization.
Actuating macrophage polarization for acceleration of wound healing must be done in a narrow time interval, beginning from the peak time of M1. Too early and too strong treatment of this type may slow down wound cleaning and lead to chronic inflammation. Delayed treatment may have too small of an effect. For the particular parameter values chosen here, the optimal actuating time is between 0.7 and 3.1 days. The shortest observed wound healing time was 11.37 days for a ∼15.5% reduction in wound healing time. This is comparable to the results in Liang et al. [53], where time for wound closure was reduced by ∼12-18%.
To our knowledge, this is the first work in search of an optimal spatio-temporal regime of wound treatment, which can be tested experimentally. We also note that we are only presenting one method of intervention. In the future, this study can be expanded to include additional intervention strategies targeting other biological processes. Thus, this naive modeling approach may help to predict optimal regimes for various treatments with known distinct actions. Combining approaches could potentially lead to unprecedented reductions in wound healing times. We believe this is the first step toward designing smart treatment and development of algorithms for smart wound healing devices.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
KZ provided context and background, significantly contributed to the model, performed numerical analysis, prepared images, and wrote parts of the original draft. JX performed numerical analysis, prepared images, and parts of the original draft. MG contributed to conceptualization, supervised the work, reviewed, and edited the article. All authors contributed to the article and approved the submitted version.