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–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–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–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 (pro-inflammatory) 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–27]. There are several substances that can induce macrophage polarization from the pro-inflammatory to the anti-inflammatory subtype [28–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.

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₀ make the treatment less effective: wound healing time increases as t₀ increases. In our simulations for t₀ > 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 time-window 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, 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₀. 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 5 shows how wound healing time depends on the duration of treatment. The longer the treatment time Δt is, the more effective it is at accelerating wound healing. However, for Δt larger than 2–3 days, wound healing time approaches a lower bound. This means that further prolongation of the treatment has minor effects. The plots for other values of r_p are shown in Supplementary 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₀ 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

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:

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₀ (left of the first dotted line) lead to too slow wound cleaning. Regimes with very large t₀ (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₀, 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₀,Δ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₀ = 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.

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 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₀, Δ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₀ 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.

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.

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

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.

Research was sponsored by the Office of Naval Research and the DARPA Biotechnologies Office (DARPA/BTO) and was accomplished under Cooperative Agreement Number DC20AC00003.

The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Office of Naval Research and the DARPA Biotechnologies Office (DARPA/BTO) or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.