Traveling wave of inflammatory response to regulate the expansion or shrinkage of skin erythema

Many skin diseases show circular red lesions on the skin, called erythema. Erythema is characterized by the expansion of its circular area solely from local stimulation. A pathological inflammatory response caused by the stimulation persistently increases inflammatory mediators in the dermis, whereas a normal inflammatory response transiently increases mediators, resulting in the shrinkage of the erythema. Although the diffusion of mediators theoretically reproduces the expansion, how the inflammatory response expands or shrinks the erythema remains unknown. A possibility is positive feedback, which affects mediator production and can generate two distinct stable states (i.e., inflamed and noninflamed), referred to as bistability. Bistability causes a state transition either from the noninflamed to inflamed state or the reverse direction by suprathreshold stimulation. Additionally, the diffusion selectively causes state transition in either direction, resulting in spatial spread of the transited state, known as the traveling wave. Therefore, we hypothesize that the traveling wave of the inflammatory response can account for both the expansion and shrinkage. Using a reaction-diffusion model with bistability, we theoretically show a possible mechanism in which the circular inflamed area expands via the traveling wave from the noninflamed to the inflamed state. During the expansion, the boundary between the inflamed and noninflamed areas moves at a constant velocity while maintaining its concentration gradient. Moreover, when the positive feedback is weak, the traveling wave selectively occurs from the inflamed to noninflamed state, shrinking the inflamed area. Whether the inflamed area expands or shrinks is mainly controlled by the balance of mediator concentration between the noninflamed and inflamed states, relative to the threshold. The traveling wave of the inflammatory response provides an experimentally testable framework for erythema expansion and shrinkage, thereby contributing to the development of effective treatments, including probiotics.

Erythema is caused by various pathogenic factors, such as physical stimulation, chemical drugs, and bacterial infections (Fig 1Ai) [1]. These factors lead to inflammatory response, that is, the secretion of inflammatory mediators, such as cytokines (e.g., TNF-α and IL-1β) and histamine (Fig 1Aii) [3][4][5]. These mediators increase blood volume in the dermal blood vessels (Fig 1Aiii) [6,7]. The increase in the blood volume appears as erythema (Fig 1Aiv), which typically appears over a few millimeters solely by local transient stimulation (e.g., in a few minutes) [8] and expands to a few centimeters in a few days [1,8]. During expansion, the lesions are well-circumscribed (i.e., thick red-colored) in some diseases or skin conditions ( Fig 1B) and poorly circumscribed (i.e., light red-colored) in others (Fig 1C) [9]. The expansion of multiple erythemas leads to their fusion [1]. Autonomous expansion is an indispensable event during disease progression. In a pathological case, inflammatory response persists, and concentrations of mediators fail to return to the original levels [10][11][12][13]. In contrast, the inflammatory response in the healthy skin (normal inflammatory response) initiates a temporal increase in the level of mediators, which returns to original levels [10][11][12][13]. Controlling the inflammatory response to suppress the expansion velocity and further shrink the erythema can offer indications for developing effective treatments. Nevertheless, how the inflammatory response controls the expansion and shrinkage of erythema remains unclear because of the experimental difficulty in detecting spatiotemporal dynamics of inflammatory mediators in the dermis.
In recent years, mathematical models and computer simulations have predicted the spatiotemporal dynamics in the dermis [14][15][16][17][18][19]. One model incorporated experimentally known biochemical or transcriptional regulation of mediators and their intercellular diffusion and highlighted the expansion of a well-circumscribed lesion [14]. This type of model incorporates the reaction and diffusion of molecules, referred to as the reaction-diffusion model [20,21]. Another reaction-diffusion model incorporating self-activation, that is, positive feedback of histamine, has also shown expansion [15]. These models show that the diffusion of mediators could cause erythema expansion. The common regulation of mediators responsible for the expansion in these two models is the positive feedback. Experiments support the positive feedback in the dermis, that is, the inflammatory mediators activate NF-κB signaling, and their production is induced in response to NF-κB activation [22]. Mathematically, positive feedback can generate two distinct and stable states (i.e., inflamed and noninflamed states) called bistability. Bistability causes a persistent transition between the noninflamed and inflamed states by a suprathreshold stimulation [13]. Diffusion and bistability selectively cause transition from one (e.g., noninflamed) state to another (e.g., inflamed), resulting in the spatial spread of the state transitions, referred to as the traveling wave [20,23]. Furthermore, weak positive feedback selectively causes a reverse transition, e.g., from inflamed state to noninflamed state, resulting in a traveling wave in the reverse direction [20,23]. Therefore, we hypothesized that the diffusion and bistability of inflammatory mediators could account for both expansion and shrinkage.
This study develops a bistable reaction-diffusion model to determine whether and how diffusion and bistability can cause expansion and shrinkage. The expansion of the inflamed area appears as the traveling wave of the transition from the noninflamed to inflamed state. We further demonstrate that diffusion and bistability can shrink the inflamed area through a traveling wave of a reverse transition from the inflamed to noninflamed state, depending on the strength of the positive feedback.

Methods
To analyze the expansion of erythema, we formulated a reaction-diffusion model [20,21,23]. Because mediator-secreting cells (e.g., immune cells and keratinocytes) do not exhibit spatial localization in the dermis [1,4,24], we assumed a homogeneous distribution of these cells in the two-dimensional space of the dermis along the skin surface ( Fig 1A). We first formulated the following equation by introducing the observed biochemical or transcriptional regulation of the inflammatory mediator's concentration (p) in the intracellular and extracellular  (3) represented by the production rate of mediators ( dq dt ) as a function of the concentration (q). Two filled circles represent stable steady states (S NI , S I ), whereas the hollow circle indicates an unstable steady state (S T ). This system can switch between the stable states of low (S NI )-and high (S I )-concentration depending on the perturbation, such as initial stimulation or diffused mediators (a = 2.14, b = 0.05).
https://doi.org/10.1371/journal.pone.0263049.g001 environments into an ordinary differential equation: The first, second, and third terms represent the induction of own production (i.e., the positive feedback) [22], basal secretion [25], and degradation [26], respectively. Here, α, n, K M , β, and γ denote the maximum production rate, Hill coefficient of the cooperativity, threshold of production, basal secretion rate, and degradation rate, respectively. The values of these parameters can depend on the skin condition. For example, experiments have suggested that the maximum production rate (α) of one type of mediator, IL-1β, increased with the deterioration of skin microbiome [27], and that the basal secretion rate (β) of IL-1β increased with the deficiency of the skin barrier integrity [25].
Then, we introduced the diffusion to formulate a reaction-diffusion equation: where D and Δ in the fourth term denote the diffusion coefficient and the Laplacian operator , respectively [28]. Eq (2) becomes identical with the previous model incorporating both inflammatory mediator and its substrate as variables [14] when the substrate is assumed to be in a steady state (See Appendix A2 in S1 Appendix for a detailed derivation). We set n (Eq 2) to 2 to introduce the simplest form of the cooperativity required for the bistability.
Because there is no quantitative information on the other kinetic parameter values (α, β, γ, and K M ) and diffusion coefficient (D), we investigated the model dynamics for a wide range of parameters, such that the model exhibits bistability. For this purpose, we non-dimensionalized Eq (2) by normalizing the variables and parameters as follows (See Appendix A1 in S1 Appendix for a detailed derivation): where q, t, a, b, and d are the normalized concentration of mediator, normalized time, normalized maximum production rate, normalized basal secretion rate, and normalized diffusion coefficient, respectively. We analytically determined the range of a and b such that the model exhibits bistability in the absence of diffusion ( Fig 1D); there are two stable states given by low and high concentrations, corresponding to the noninflamed (S NI in Fig 1E) and inflamed (S I in Fig 1E) state, respectively. We assumed that an area with a concentration of S I in the dermis appears as erythema on the skin surface. In this setting, the model has one unstable steady state, corresponding to a threshold concentration (S T in Fig 1E) for the transition between the two stable states. S NI is maintained for a subthreshold perturbation (i.e., below the concentration of S T ), whereas it transits to S I for a suprathreshold perturbation (i.e., above the concentration of S T ). S I also transits to S NI when the concentration decreases below S T .
Finally, as an initial condition of the model simulation (Eq 3), we referred to the physiological condition at the onset of erythema, where one or a few small (~1mm) inflamed areas exhibited a concentration of mediators above the threshold (S T ) [1]. In contrast, the surrounding areas exhibited a concentration of mediators below the threshold in the dermis [1]. Based on these observations, for each inflamed area, we set a circular area of q > S T and the surrounding area of q < S T , which are given by the two-dimensional Gaussian distribution (Fig 2A, t = 0). Given this initial condition, numerical simulation of Eq (3) was performed in two-dimensional  where Δt, Δx, and Δy were chosen to satisfy Von Neuman stability. We confirmed that the obtained results were barely influenced by the choice of the temporal discretion size Δt ( Fig  2D). A simulation code written in C language is available from GitHub: https://github.com/ MakiSudo/Travelingwave_Simulation/blob/bc2c10ddd5eff8db374b0804e11a63ef3c0e766a/ Simulationcode.c.

Diffusion and bistability can cause expansion of circular inflamed area
We examined whether diffusion and bistability can cause expansion of the erythema in the model. The model simulations showed that a circular inflamed area was initially caused by a transient and local perturbation to the mediator's concentration and subsequently expanded centrifugally over time (Fig 2A), consistently with the expansion of erythema (Fig 1B). During the expansion, the inflamed area maintained a steep gradient of concentration at the boundary ( Fig 2B) and increased the diameter at a constant rate (velocity) over time (Fig 2C and 2D). The mediator level was persistently high (S1A and S1B Fig in S1 Appendix), consistently with the pathological inflammatory response [13]. This model further reproduced the fusion of multiple inflamed areas (Fig 2E). Even when the ratio of the inflamed and noninflamed concentration (S I /S NI ) was smaller according to a decrease in the maximum production rate (a), the inflamed area similarly expanded ( Fig 2F) with a steep boundary gradient ( Fig 2G) at a constant velocity (Fig 2H and 2I) and fused (Fig 2J). We then examined whether diffusion is necessary for expansion. Without diffusion (i.e., d = 0), an inflamed area appeared; however, this area did not expand and remained constant over time (S1C Fig in S1 Appendix). We next examined whether bistability is necessary for expansion. When bistability was lost by further decreasing a, a local transient stimulation caused an inflamed area, but the area did not expand (S1D Fig in S1 Appendix). Similar results were obtained when the bistability was lost by decreasing the basal secretion rate (b) (S1E Fig in S1 Appendix). Thus, these results confirmed that diffusion and bistability could cause the expansion of the circular inflamed area.

Expansion caused by a traveling wave of the transition from a noninflamed to inflamed state
The expansion follows spatiotemporal changes in inflammatory mediator concentration. First, an initial suprathreshold perturbation locally induces an inflamed area (Fig 2A, 2B, 2F and 2G, time = 0). In this area, the production of mediators increases by positive feedback. The produced mediators diffuse to the adjacent noninflamed area. In the noninflamed area, the diffused mediators are large enough to become a suprathreshold perturbation, causing a selective transition from the noninflamed state (S NI ) to the inflamed state (S I ) at the boundary between the inflamed area and the noninflamed area (Fig 2A, 2B, 2F and 2G, e.g., time = 150). This series of events, that is, positive feedback of production, diffusion, and state transition in the adjacent area, occurs in each position and propagates to the surrounding noninflamed area. Thus, the inflamed area expands as the traveling wave, while maintaining the velocity (Fig 2D  and 2I) and a gradient at the boundary (Fig 2B and 2G). Therefore, diffusion and bistability cause the traveling wave of selective transition from the noninflamed to the inflamed state, resulting in the expansion.

Control of bistability or diffusion to suppress the expansion velocity
We examined whether the expansion velocity of the inflamed area can be suppressed by controlling diffusion and bistability. The expansion velocity monotonically decreased with a decrease in the diffusion coefficient (d) (Fig 3A). It also decreased with a decrease in the maximum production rate (a) or basal secretion rate (b) controlling the positive feedback activity (Fig 3B-3E and S2A-S2D Fig in S1 Appendix). Unlike the dependence on d, the velocity continuously decreased and fell below zero at a threshold value of a and b (Fig 3B, 3C and 3F). Note that the parameter value did not affect the gradient (Fig 3G-3I). Thus, these results show that the expansion velocity was suppressed by diffusion and positive feedback activity for a given bistability.

Erythema shrinkage by controlling bistability
Below the threshold value of a or b, the expansion velocity became negative, where transition selectively occurred from the inflamed (S I ) to the noninflamed (S NI ) state ( Fig 3B, 3C and 3F). This traveling wave resulted in the shrinkage of the inflamed area (Fig 3B, 3C, 3J and S2E-S2G Fig in S1 Appendix). During the shrinkage, the mediator level at the initial stimulation transiently increased to the inflamed state and then returned to the original noninflamed state (S2H Fig in S1 Appendix), consistently with the normal inflammatory response [13]. Therefore, these results showed that diffusion and bistability accounted for the normal inflammatory response leading to the shrinkage as well as the pathological response leading to the expansion.

Balance between the inflamed and noninflamed state concentrations determines expansion or shrinkage
Finally, we clarify how bistability controls the expansion and shrinkage (Fig 4A). To theoretically formulate the velocity of the traveling wave using the model parameters [20], we approximated Eq (3) to under an assumption that 1 S T n þ1 is approximated to a constant A (See Appendix A3 in S1 Appendix). Mathematically, the velocity is approximately determined by the diffusion coefficient (d) and the concentrations at the three steady states (i.e., S NI , S T , The theoretical velocity is proportional to the square root of the diffusion coefficient (d), and more interestingly, to the difference between the concentrations of the stable (S NI and S I ) relative to unstable (S T ) states. For the range of parameter values showing bistability, the theoretical velocity agreed with the simulated velocity in the sign (Figs 3F and 4A) and the approximate value, except for the parameter values at the velocity of zero (S3 Fig in S1 Appendix). When S T is closer to S NI than to S I , the velocity is positive, indicating the expansion of the inflamed area (Fig 4A and 4B). When S T is at an equal distance from S NI and S I , given a decrease in the maximum production rate (a), the velocity is suppressed toward zero (Fig 4A  Fig 3. Expansion velocity is controlled by diffusion and positive feedback for given bistability. (A-C and 4C). The velocity is negative, indicating the shrinkage of the inflamed area when S T is closer to S I than to S NI (Fig 4A and 4D). Similar results were obtained by decreasing the basal secretion rate (b). Therefore, depending on whether the threshold concentration (S T ) is closer to the noninflamed state (S NI ) or inflamed state (S I ), erythema expands or shrinks, respectively

PLOS ONE
( Fig 5). The balance of mediator concentrations (i.e., how close the threshold is to the noninflamed or inflamed state) further determines the velocity of expansion or shrinkage.

Diffusion and bistability cause both expansion and shrinkage of erythema as the traveling wave
Erythema is characterized by the expansion of its circular area. In pathological skin showing expansion, inflammatory response persists and the concentration of mediators fails to return to the original level [10][11][12][13]. Recent theoretical studies have shown two independent mechanisms of how the pathological inflammatory response causes expansion [14,15]. One is the diffusion of inflammatory mediators and the other is positive feedback, which can generate bistability. However, the mechanism of how diffusion and bistability can cause expansion remains unknown. Furthermore, how we can reduce the expansion velocity and further shrink the erythema has not yet been identified. In this study, we theoretically show that diffusion and bistability can synergistically cause not only expansion (Fig 2) with the pathological inflammatory response (S1B Fig in S1 Appendix) but also shrinkage (Fig 3; S2E-S2G Fig in S1 Appendix) with the normal inflammatory response (S2H Fig in S1 Appendix) by the traveling wave. Whether the inflamed area expands or shrinks is determined by the concentration balance between the noninflamed (S NI ) and inflamed (S I ) state relative to the threshold (S T ; Eq 5). Expansion occurs when S T is closer to S NI than to S I , whereas shrinkage occurs when S T is closer to S I than to S NI in a wide range of parameters controlling the inflammatory response, including positive feedback, basal secretion, and degradation (Fig 4; S3 Fig in S1 Appendix). An interesting future study would be the analysis of whether the balance captures such inflammatory wave dynamics in other bistable systems, with more complex biochemical reactions [12,14,15]. Therefore, the balance of bistable states could provide an experimentally testable framework for the normal and pathological inflammatory responses (Fig 5).

Expansion of a well-and poorly-circumscribed erythema
Erythema expands with well-circumscribed lesions ( Fig 1B) or with poorly circumscribed lesions ( Fig 1C) depending on diseases and skin conditions [9]. A well-circumscribed lesion indicates an inflamed area clearly distinguished from the surrounding noninflamed area [29], which may appear in situations such as a sharp gradient at the boundary and/or a large ratio of the inflamed and noninflamed concentration (S I /S NI ) in the dermis. In contrast, a poorly circumscribed erythema, which is difficult to distinguish from the surrounding noninflamed area [29], may appear in situations such as a shallow gradient or a small concentration ratio (S I /S NI ) in the dermis. Some of the expansion of inflamed areas in our model appear with the steep gradient and a large ratio of the inflamed and noninflamed concentration (S I /S NI ) in the dermis (Fig 2A and 2E), accounting for the possible situations of well-circumscribed lesions. Other expansions appear with a small ratio (S I /S NI ) in the dermis (Fig 2F and 2J), consistently with the one possible situation of poorly circumscribed lesions. The consistency can be experimentally verified by checking whether a poorly circumscribed erythema has a small concentration ratio and expands at a constant velocity (e.g., Fig 2I). Alternatively, a well-circumscribed lesion in a deep layer of the dermis may appear to be a poorly circumscribed lesion on the skin surface [29]. Our framework for the expansion of both well-and poorly-circumscribed lesions provides a new perspective on how inflammation in the dermis appears as erythema on the skin surface.

Possible relevance to biological factors in the skin
The expansion velocity was suppressed by decreasing the model parameters (Fig 3), which can depend on the skin condition. Experiments have shown that the maximum production rate (α in Eq 2, a in Eq 3) and basal secretion rate (β in Eq 2, b in Eq 3) are lower in healthy skin than in pathological skin with a deterioration of the skin microbiome [27] and in skin with deficiency of the skin barrier integrity [25], respectively. Measuring the expansion velocity under different skin conditions will reveal the relation between the model parameters and the skin conditions, thereby potentially providing possible treatments to lower the maximum production rate or the basal secretion rate. For example, probiotics that improve the skin microbiome composition significantly lower the maximum production rate [30]. Additionally, probiotics improve the skin barrier integrity [30], which is expected to lower the basal secretion rate. Thus, probiotics can lower the maximum production rate and the basal secretion rate, thereby possibly reducing the expansion velocity and shrinking erythema. Further study of the relationship between the expansion velocity and the skin condition will offer further insights helpful in developing more effective treatments of erythema.

Limitation of the present model
In pathological inflammatory response, erythema expands to a certain size and eventually autonomously disappears [1]. Unlike shrinkage in normal inflammatory response (Fig 3J;  S2E-S2H Fig in S1 Appendix), the disappearance often shows a decrease in intensity (i.e., redness) without changing the diameter of erythema [14]. Such disappearance could not be reproduced by the present model of inflammatory mediator alone; it may require other factors. Some of the responsible factors are anti-inflammatory mediators, such as IL-10 and TGF-β [5,6]. The anti-inflammatory mediators are produced by inflammatory mediators during the development of erythema, and inhibit the production of inflammatory mediators [22]. A previous mathematical model incorporating the interaction of these mediators accounted for the temporal evolution (i.e., decrease) of the inflammatory response [13], but not the spatiotemporal evolution (i.e., autonomous disappearance). Thus, future studies should extend our model to incorporate the interaction and examine how anti-inflammatory and inflammatory mediators synergistically control the disappearance of erythema.

Conclusions
In this paper, we demonstrated that diffusion and bistability could cause expansion and shrinkage as a traveling wave. Furthermore, the positive feedback activity of mediator production regulates the transition between the noninflamed and inflamed states, thereby determining whether the inflamed area expands or shrinks (Fig 5). Moreover, regulating the balance of mediator concentration between noninflamed and inflamed states provides an experimentally testable framework for the spatiotemporal evolution of erythema, which can help in the development of effective treatments.