Stability Analysis of a Model for Foreign Body Fibrotic Reactions

Abstract

Implanted medical devices often trigger immunological and inflammatory reactions from surrounding tissues. The foreign body-mediated tissue responses may result in varying degrees of fibrotic tissue formation. There is an intensive research interest in the area of wound healing modeling, and quantitative methods are proposed to systematically study the behavior of this complex system of multiple cells, proteins, and enzymes. This paper introduces a kinetics-based model for analyzing reactions of various cells/proteins and biochemical processes as well as their transient behavior during the implant healing in 2-dimensional space. In particular, we provide a detailed modeling study of different roles of macrophages (𝑀Φ) and their effects on fibrotic reactions. The main mathematical result indicates that the stability of the inflamed steady state depends primarily on the reaction dynamics of the system. However, if the said equilibrium is unstable by its reaction-only system, the spatial diffusion and chemotactic effects can help to stabilize when the model is dominated by classical and regulatory macrophages over the inflammatory macrophages. The mathematical proof and counter examples are given for these conclusions.

1. Introduction

Recently, intensive research efforts have been focusing on developing mechanistic computational models for wound healing related processes. Wound healing is a very complicated biochemical and biophysical phenomenon, with many facets and subprocesses, including the inflammatory response process, angiogenesis as well associated fibrotic reactions. Many cells, enzyme, growth factors, and proteins participate at different stages of the wound healing reactions, and they form a network of signaling pathways that in turn leads to inflammatory, angiogenesis, and fibrotic reactions. We refer to the review by Diegelmann and Diegelmann and Evans 2004 [1] for a brief review of the recent scientific work.

As a subarea of general wound healing research, healing processes involved in medical implantations are of significant application for modern medicine [2–4]. It is commonly accepted that implants may cause foreign body reactions that are initiated with implant-mediated fibrin clot formation, followed by acute inflammatory responses [4, 5]. The inflammatory chemokines released by adherent immune cells serve as strong signals for triggering the migration of macrophages and fibroblasts from the surrounding tissues and circulation toward the implant surface [5]. The implant-recruited fibroblasts consequently synthesize chains of amino acids called procollagen, a process that is activated by growth factors, including in particular type-𝛽 transforming growth factor (TGF-𝛽) [6, 7] to become collagen, the dominant ingredient of the extracellular matrix (ECM) [8]. These processes may, however, differ slightly between dermal wound healing and implantation when it comes to specific activation and inhibition loops of reactions.

Among inflammatory cells, macrophages (𝑀Φ) are found to reside in the wound [9]. The roles of macrophages are multiple and stand prominent in the activations and inflammations during implantation. 𝑀Φ are known to remove damaged tissue and foreign debris via phagocytosis. In addition, 𝑀Φ often release a variety of chemokines to recruit other cell types, such as fibroblasts, which participate in the remodeling of ECM. The specific roles of 𝑀Φ vary significantly at different stages of healing process. Work by Mosser and Edwards 2008 [10] has shown there to be at least 3 phenotypes of 𝑀Φ, each of which displays a different functionality. Classically activated 𝑀Φ represent the effector 𝑀Φ that are produced during cell-mediated immune responses. Two signals, interferon-𝛾 and tumor-necrosis factor-𝛼, give rise to these effector 𝑀Φ which have enhanced microbicidal or tumoricidal capacity and secrete high levels of proinflammatory cytokines and mediators. Assisted in part by the production of transforming growth factor type 𝛽 (TGF-𝛽), the clearance of apoptotic inflammatory, as well as noninflammatory cells by classical 𝑀Φ, can lead to an inhibition of inflammation [11, 12]. Wound healing 𝑀Φ (or inflammatory 𝑀Φ) can develop in response to innate or adaptive signals through interleukin-4. In turn, interleukin-4 stimulates arginase activity in 𝑀Φ, allowing them to convert arginine to ornithine, a precursor of polyamine and collagen, thereby contributing to the production of extracellular matrix (ECM) [13]. Regulatory 𝑀Φ can also arise during the later stages of adaptive immune responses, the primary role of which dampen, the immune response and limits inflammation through production of interleukin-10 [14]. Although all three phenotypes were observed experimentally within the dermal wound healing context, the phagocyte biomaterial interactions are known to be similar here for foreign body reactions.

While experiments are still the main stay in the studying of wound healing related process, significant progress has also been made in detail predictive modeling based on biochemical and biophysics principles. For dermal wound healing, basic reactions were first considered in studies by Dale et al. 1996 [15], 1997 [16]; Dallon et al. 2001 [17] and many others. Their models incorporated the key features of kinetics which are essential to dermal wound healing. Their results have successfully described the dynamics and compare favorably with experiments, in terms of healed ECM fiber ratio, spatial orientation, and other features. Recently, the work of Schugart et al. 2008 [18] and Xue et al. 2009 [19] further included angiogenesis equations to the healing process and examined the positive effects of increased oxygen level in accelerating the healing and closure of open wound, suggesting new insights for the healing. Furthermore, a wound healing model based more on cell migration was considered in Arciero et al. 2011 [20].

Inflammatory reactions are important to wound healing as they activate many key agents for the healing process, however, prolonged inflammation may cause excessive scars and chronic wounds. Through interactions between immune mediators, phagocytes in the blood and tissue, the acute inflammatory response was modeled and analyzed by reduced compartmental models in Reynolds et al. 2006 [21] and Day et al. 2006 [22]. Closely related to dermal wound healing and implantation, atherogenesis in blood vessels was modeled by continuum equations in Ibragimov et al. [23]. The concept of debris and phagocytosis in [21, 22, 24] is analogous to our current model, which assumes that the digestion of dead cells (or tissues) initiates the entire healing process. Further, addition of stem cells can create a new dimension to the healing and implantation process; we mention Lemon et al. 2009 [24] for their new mathematical tool in providing quantitative analysis for this growing field.

Our primary goal in this paper is to use computational modeling to study the fibrotic reaction process following implantation with specific attention given to the effects caused by varying the mix of different phenotypes of 𝑀Φ. Our modeling results indicate trends for these variations, serving as a plausible clue for developing new experiments.

The main mathematical contribution of this paper is as follows. The nonzero equilibrium of our model represents an inflamed state. If it is linearly stable in terms of the corresponding ODE system (the reactions network of the model), then it is also stable for the full system (which includes spatial diffusion and chemotaxis). In other words, spatial effects cannot destabilize the equilibrium if it is stable in its pure reactions. However, even if the equilibrium is unstable by its reaction system, the spatial diffusion and chemotactic effects can help to stabilize the equilibrium under several conditions. These conditions suggest the need for the model to be dominated by classical and regulatory macrophages over the inflammatory macrophages. The mathematical proof and counter examples are given for these results.

We organize the paper as follows. In Section 2, we introduce the model and the modeling considerations. In Section 3, we discuss the spatially uniform equilibria and their stability in relation to the ODE system representing the reaction system without spatial variations. In Section 4, we prove that if the equilibrium is stable in ODE sense, then it is stable for the full system with respect to any small spatial perturbation measured in 𝐿². In Section 5, we provide a set of sufficient conditions for the equilibrium to be stable for the full system, which allows us to explicitly give a counter example where PDE solutions can be conditionally stable without requiring stability in the ODE system. A brief summary and discussion are presented in Section 6.

2. Modeling Based on Chemical Kinetics Equations

Our foreign body reaction model is partially from the mass-action kinetics framework developed by Schugart et al. 2008 [18], which modeled wound healing under oxygen pressure. In medical implantation processes, new kinetics of 𝑀Φ reactions were added in the framework. The main biological question that we hope to address is the variance between tissue responses at different percentages of classical, inflammatory, and regulatory 𝑀Φ cells during foreign body fibrotic reaction processes. We model the following: 𝜕𝐷𝜕𝑡=𝐷^𝑑𝜕𝐶∇²𝐷−𝑓₀𝜆₁𝑀𝐷+𝑓₀𝜆₃𝑀,(1)𝜕𝑡=𝐷^𝑐𝜕𝐹∇²𝐶+𝑓₁𝐷+𝑓₂𝜆₃𝑀−𝑓₃𝜆₂𝑀𝐶−𝑓₄𝐶,(2)𝜕𝑡=𝐷^𝑓+𝑎₂∇²𝐹−𝜒₀∇⋅(𝐹∇𝐶)+𝑎₁𝜆₁𝑀𝐹𝑎²−𝑎₃𝑎₂−𝐹𝐹₀+𝑎¹²𝜕𝑀𝐶𝐹𝐻𝐹₀−𝐹,(3)𝜕𝑡=𝐷^𝑚−𝑎₀∇²𝑀−𝜒₁∇⋅𝑀𝐻𝑀₀−𝑀∇𝐶𝜕𝐸𝑀+𝑎₁₁𝐶𝑀𝐻𝑀₀−𝑀,(4)𝜕𝑡=𝐷^𝑒∇²𝐸−∇⋅Φ+𝑎₁₆𝐹1−𝐸𝐸₀,(5) where ∇²=∇⋅∇, and the vector field Φ=𝐵𝐷^𝑓𝐹₀𝐸∇𝐹+𝐵𝜒^𝑗𝐹₀𝐸𝐹𝐻𝐹₀−𝐹∇𝐶,(6) and all coefficients are positive. The form of the logistic terms in (3) is for representing biological meanings of the coefficients.

In the system (1)–(5), the debris cell population 𝐷 represents dead tissue cells following implantation. Abnormal white blood cells and molecules caused by the surgery are also included in this debris term, which is assumed to be the initiation point of reactions. We assume that they are digested by 𝑀₁-classical 𝑀Φ and that 𝑀₃-inflammatory 𝑀Φ contribute to the accumulation of debris during the healing process (as modeled in (1)).

The chemoattractant consists mainly of various forms of growth factors including tissue growth factors type 𝛽 (TGF𝛽) released during the tissue injury. The chemoattractant field 𝐶 is enhanced by the presence of debris and 𝑀₃-inflammatory 𝑀Φ cell, but is inhibited by 𝑀₂-regulatory 𝑀Φ cells. In (2), we assume that cell spatial migration occurs through diffusion and chemotactic migration based on the gradient field of 𝐶.

Fibroblast density 𝐹 represents a main cell type in secreting collagen (a major component of ECM). Fibroblast proliferation and collagen synthesis are upregulated by the chemoattractant gradient field 𝐶. Thus fibroblast population 𝐹 (shown in (3)) can be approximated by a chemically enhanced logistic growth 𝐹(1−(𝐹/𝐹₀)) with a threshold 𝐹₀, along with its diffusion in space modeled by 𝐷_𝑓∇²𝐹, chemotactic migration by −𝜒₀∇⋅(𝐹∇𝐶) and its natural decay according to time as shown in (3). New experimental data also shows autocrine upregulation of fibroblast by TGF𝛽 without chemotaxis [25]; this effect is also included in the modeling. The term 𝑎₃𝐹 is the decaying factor.

Macrophage density, 𝑀, is the summation of 𝑀₁-classical 𝑀Φ, 𝑀₂-regulatory 𝑀Φ, and 𝑀₃-inflammatory 𝑀Φ. We assume that they each take on a proportion 𝜆₁, 𝜆₂, and 𝜆₃ of 𝑀Φ, respectively. Each phenotype 𝑀_𝑗, 𝑗=1,2,3, may take a different share of 𝑀Φ at different stages of foreign body fibrotic reactions. However, our model simplifies the situation in that (a) the proportions 𝜆₁, 𝜆₂, and 𝜆₃ for different phenotypes of 𝑀Φ are fixed, and (b) the total 𝑀Φ population is set to share one common biochemical reaction equation (4), since its basic biochemical properties are similar. The proliferation of 𝑀Φ at the field is through diffusion and migration upregulated by the chemotactic gradient field 𝐶, but the production does reach a limiting value once the 𝑀Φ population reaches its saturation of 𝑀₀. 𝑀Φ cell apoptosis and proliferation caused by the direct interaction with chemoattractants are also assumed.

Finally in (5), fibroblasts secrete procollagen which is then activated by the chemoattractant TGF𝛽s into collagen (or ECM) represented by the quantity 𝐸. We also incorporate the effects of ECM diffusion, fibroblast movement, chemotactic migration, and ECM saturation in mass-action law. In all discussions, 𝐻 is the Heaviside function, and 𝑀₀ is the 𝑀Φ saturation level.

We assume in our implant model that the computational domain is large enough and also the cell changes are slow enough (measured in days) that there is no significant boundary flux, allowing us to take homogeneous Neumann boundary conditions as a reasonable approximation.

Definition 1. Let us define inflammatory equilibrium as a strictly nonzero constant vector 𝑈_𝑒 in 5-dimensional space 𝑈_𝑒=(𝑑_𝑒,𝑐_𝑒,𝑓_𝑒,𝑚_𝑒,𝑒_𝑒) with 𝑑_𝑒>0,𝑐_𝑒>0,𝐹_𝑜≥𝑓_𝑒>0,𝑀_𝑜≥𝑚_𝑒>0,𝑒_𝑒=𝐸₀>0, which solves system of the equations (1)–(5).

Remark 2. In the case of a no-flux boundary condition, the spatially uniform steady state is often used when modeling inflammatory response in tissue (see e.g., [23, 26] and reference therein). A physically realistic, nonnegative set of equilibriums can easily be obtained by letting the RHS of the original system (1)–(5) equal to zero. It is natural to define the trivial (zero) equilibrium as ground or healthy state and study its stability. Instability of the ground state is usually interpreted as unfavorable development of the disease. In this paper we take a different approach and are interested in analyzing the stability of the abnormal/inflammatory equilibrium which is nonzero for all five components of the unknown. This equilibrium can be stable or unstable depending on the parameters of the model. In this case, instability of the equilibrium does not necessarily mean an unhealthy response of the immune system. An instability of a nonzero equilibrium can lead to a ground healthy state (best case scenario), to another steady state (uncertain developments), or to infinity (acute development). If in contrary, the perturbation of 𝑈_𝑒 is linearly stable and vanishes at time infinity, then 𝑈_𝑒 can be interpreted as sustainable. All these make linear stability analysis very appealing from both a theoretical and applied point of view. It is worth mentioning that from a biological point of view, a strictly positive steady state 𝑈_𝑒 can be transitioned from some other nonstrictly positive state. We believe that this type of interpretation of the inflammatory equilibrium stability conditions is logical and presents an example of a sustainable wound which does not heal over the course of a long time period (see [19–21]). An indirect analogy of such an inflammatory (chronically) stable equilibrium has been introduced and applied for studying biological dynamic system in virology for some years (see e.g., [27]). At this stage of the research, we are studying stability of the strictly positive state 𝑈_𝑒 mostly as a model of inflammatory equilibrium, without analysis of its genesis. As commonly occurs in biomedical research, the mathematical model can often provide nonintuitive insights into dynamics of inflammatory responses in the wound healing processes and can suggest new avenues for experimentation. In the forthcoming sections, sufficient conditions on the parameters of the system of the equation guarantee stability of nonzero equilibrium.

2.1. Linearized System

Let perturbation near this equilibrium be as following: 𝑑=𝐷𝑚−=𝑑_𝑒𝑀,−𝑐=𝑚_𝑒𝐶−𝑐_𝑒,𝑓=𝐹−𝑓_𝑒,,𝑒=𝐸−𝑒_𝑒.(7) Denote vector field of the perturbation by 𝑣(𝑥,𝑡)=(𝑑,𝑐,𝑓,𝑚,𝑒). Then the linearized system for 𝑣(𝑥,𝑡) will take the following form: 𝜕𝑑𝜕𝑡=𝐷^𝑑𝜕𝑐∇²𝑑−𝑏₁₁𝑑−𝑏₁₄𝑚𝜕𝑡=𝐷^𝑐𝜕𝑓∇²𝑐−𝑏₂₁𝑑−𝑏₂₂𝑐−𝑏₂₄𝑚,𝜕𝑡=𝐷^𝑓𝜕𝑚∇²𝑓−𝜒_𝑓∇²𝑐−𝑏₃₂𝑐−𝑏₃₃𝑓−𝑏₃₄𝑚,𝜕𝑡=𝐷^𝑚𝜕𝑒∇²𝑚−𝜒_𝑚∇²𝑐−𝑏₄₂𝑐−𝑏₄₄𝑚,𝜕𝑡=𝐷^𝑒∇²𝑒−𝜒_𝑒1∇²𝑓−𝜒_𝑒2∇²𝑐−𝑏₅₃𝑓−𝑏₅₅𝑒.(8) Here, 𝜒_𝑓=𝑓_𝑒𝜒₀,𝜒_𝑚=𝜒₁𝑚_𝑒,𝜒_𝑒1=𝐵𝐷^𝑓𝑒₀𝐹₀,𝜒_𝑒2=𝐵𝜒^𝑗𝑒₀𝑓_𝑒𝐹₀𝑏₁₁,𝑏=𝑓₀𝜆₁𝑚_𝑒,𝑏₁₂=0,𝑏₁₃=0,₁₄𝑏=−𝑓₀𝜆₃−𝑓₀𝜆₁𝑑_𝑒,𝑏¹⁵=0,₂₁𝑏₂₄=−𝑓₁,𝑏₂₂=𝑓₃𝜆₂𝑚_𝑒+𝑓₄,𝑏=−𝑓₂𝜆₃−𝑓₃𝜆₂𝑐_𝑒,𝑏₂₃=𝑏₂₅=0,₃₁𝑏₃₃=0,𝑏₃₂=−𝑎₁₂𝑓_𝑒,=−𝑎₂1−2𝐹𝑓^𝑒₀+𝑎¹²𝑐_𝑒−𝑎₃𝑏₃₄,𝑏=−𝑎₁𝜆₁,𝑏₃₅=0,₄₁𝑏=0,𝑏₄₂=−𝑎₁₁𝑚_𝑒,𝑏₄₃=0,₄₄𝑏=𝑎₀−𝑎₁₁𝑚_𝑒,𝑏₄₅=0,₅₁=0,𝑏₅₂=0,𝑏₅₃=−𝑎₁₆𝑏₅₄1−𝐸𝑒^𝑒₀,=0,𝑏₅₅=𝑎₁₆𝑓_𝑒𝐸₀.(9)

3. Spatially Uniform Equilibrium States and Linear Stability in ODE System

We now focus on equilibrium states that are uniform in space for this Neumann problem. By removing the spatial variations, (1)–(5) reduce to the following ODE system: 𝑑𝐷𝑑𝑡=−𝑓⁰𝑑𝐶𝜆₁𝑀𝐷+𝑓₀𝜆₃𝑀,𝑑𝑡=𝑓¹𝑑𝐹𝐷+𝑓₂𝜆₃𝑀−𝑓₃𝜆₂𝑀𝐶−𝑓₄𝐶,𝑑𝑡=𝑎¹𝜆₁𝑀+𝑎₂𝐹1−𝐹𝐹₀−𝑎³𝑑𝑀𝐹+𝑎₁₂𝐶𝐹𝐻𝐹₀−𝐹,𝑑𝑡=−𝑎⁰𝑑𝐸𝑀+𝑎₁₁𝐶𝑀𝐻𝑀₀−𝑀,𝑑𝑡=𝑎¹⁶𝐹1−𝐸𝐸₀.(10)

In looking for the equilibrium of the simplified system, (10), we assume that our values are taken to be below threshold and therefore we ignore the Heaviside functions. There are several possible equilibrium states, but as it was pointed out earlier, we focus on what one can call the interior equilibrium, one in which none of the components of the equilibrium are zero. We let the right-hand side of (10) to be zero. After some algebraic work, one can obtain the following explicit formula for a unique, nonzero solution 𝑈_𝑒=(𝑑_𝑒,𝑐_𝑒,𝑒_𝑒,𝑚_𝑒,𝑓_𝑒): 𝑑_𝑒=𝑓⁰𝜆₃𝑓₀𝜆₁𝑐_𝑒,=𝑎⁰𝑎₁₁𝑒_𝑒,𝑚_𝑒=𝐸₀,=𝑓⁴𝑓₀𝜆₁𝑎₀−𝑎₁₁𝑓₁𝑓₀𝜆₃𝑓_𝑒𝑓₀𝜆₁𝑓₂𝜆₃𝑎₁₁−𝑓₃𝑎₀𝜆₂,=𝐹⁰2𝑎₂𝑎²−𝑎₃+𝑎₁₂𝑎𝑎⁰₁₁+𝐿¹.(11) Here, 𝐿₁=(𝑎₂−𝑎₃+𝑎₁₂(𝑎₀/𝑎₁₁))²+4(𝑎₂/𝐹₀)𝑎₁𝜆₁𝑚_𝑒.

Remark 3. In order for the inflammatory equilibrium to exist, it is necessary and sufficient that macrophage percentages satisfy the following: 𝑓₄𝑓₀𝜆₁𝑎₀−𝑎₁₁𝑓₁𝑓₀𝜆₃𝑓₂𝜆₃𝑎₁₁−𝑓₃𝑎₀𝜆₂>0,(12) requiring either 𝑓₂𝜆₃𝑎₁₁>𝑓₃𝑎₀𝜆₂,𝑓₄𝑓₀𝜆₁𝑎₀>𝑎₁₁𝑓₁𝑓₀𝜆₃,(13) or 𝑓₂𝜆₃𝑎₁₁<𝑓₃𝑎₀𝜆₂,𝑓₄𝑓₀𝜆₁𝑎₀<𝑎₁₁𝑓₁𝑓₀𝜆₃.(14)

Condition on the parameters in (12) says that inflammatory macrophages dominate over either regulatory or classical macrophages and are guaranteeing existence of the inflamed steady state. This point will be expounded on further in the analysis of the conditions for stability of the nonzero equilibrium state. The illustration (Figure 1) provides a visualization of the necessary macrophage phenotype parameter ranges. “Hereafter we assume that the parameters of the original model satisfy condition (12).”

Figure 1 

Illustration of the parameter range to ensure that 𝑚𝑒>0.

Turning now to satisfy the stability of the system at the equilibrium, we find the linearized system to be as follows: 𝑑𝑑𝑑𝑡=−𝑏¹¹𝑑𝑐𝑑−𝑏₁₄𝑚,𝑑𝑡=−𝑏²¹𝑑𝑓𝑑−𝑏₂₂𝑐−𝑏₂₄𝑚,𝑑𝑡=−𝑏³²𝑑𝑚𝑐−𝑏₃₃𝑓−𝑏₃₄𝑚,𝑑𝑡=−𝑏⁴²𝑑𝑒𝑐−𝑏₄₄𝑚,𝑑𝑡=−𝑏⁵³𝑓−𝑏₅₅𝑒,(15) where 𝑏₁₁=𝑏𝑓₀₂₁𝜆₁𝑚_𝑒,𝑏₁₄=−𝑓₀𝜆₃−𝑓₀𝜆₁𝑑_𝑒,𝑏₂₄=−𝑓₁,𝑏₂₂=𝑓₃𝜆₂𝑚_𝑒+𝑓₄,𝑏₃₂=−𝑓₂𝜆₃−𝑓₃𝜆₂𝑐_𝑒,=−𝑎₁₂𝑓_𝑒,𝑏₃₃=−𝑎₂1−2𝐹𝑓^𝑒₀+𝑎¹²𝑐_𝑒−𝑎₃𝑏₃₄,𝑏₄₂=−𝑎₁𝜆₁,𝑏₅₃=−𝑎₁₁𝑚_𝑒,𝑏₄₄=𝑎₀−𝑎₁₁𝑐_𝑒,=−𝑎₁₆1−𝐸𝑒^𝑒₀,𝑏⁵⁵=𝑎₁₆𝑓₀𝐸₀.(16)

Equations (32)–(39) in matrix form yields as follows: ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑚𝑓𝑑𝑐𝑒⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦^=−𝐁⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑚𝑓𝑑𝑐𝑒⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦,(17) where 𝐁 is: 𝐁=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑏¹¹𝑏₂₁00𝑏₁₄00𝑏𝑏₂₂₃₂0𝑏₂₄00𝑏₄₂𝑏₃₃𝑏₃₄000𝑏0₅₃𝑏₄₄00𝑏₅₅⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦.(18) For stability analysis, we look at the eigenvalues of matrix −𝐁; for convenience, we rearrange our equations in the following form: ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑚𝑓𝑑𝑐𝑒⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦^=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣−𝑏⁴⁴−𝑏₁₄0−𝑏₄₂00−𝑏₂₄−𝑏₁₁000−𝑏₃₄−𝑏₂₁−𝑏₂₂00000−𝑏0₃₂−−𝑏𝑏₃₃₅₃0−𝑏₅₅⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑚𝑓𝑑𝑐𝑒⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦.(19) We break −𝐁 into a 3-block and a 2-block as follows: −𝐁_𝟏=⎡⎢⎢⎢⎢⎣−𝑏⁴⁴0−𝑏₄₂−−𝑏𝑏₁₄₂₄−𝑏₁₁0−𝑏₂₁−𝑏₂₂⎤⎥⎥⎥⎥⎦,−𝐁_𝟐=⎡⎢⎢⎣−𝑏³³−𝑏₅₃0−𝑏₅₅⎤⎥⎥⎦.(20) Since det(−𝐁−𝜎𝐈)= det(−𝐁_𝟏−𝜎𝐈) det(−𝐁_𝟐−𝜎𝐈), we find the eigenvalues by looking at the eigenvalues of the 3-block, −𝐁_𝟏, and the two block, −𝐁_𝟐, separately. We also simplify by noting that with the equilibrium values found above, 𝑏₄₄=0 and 𝑏₁₄=0 s.t. det−𝐁_𝟏−𝜎𝐈=⎡⎢⎢⎢⎢⎣−−𝑏0𝜎₂₄−𝑏₁₁0−𝜎−𝑏0⁴²=−𝜎𝑏₁₁−𝑏₂₁−𝑏₂₂−𝜎⎤⎥⎥⎥⎥⎦=−𝑏₁₁+𝜎𝑏₂₂+𝜎+𝑏₂₄𝑏₄₂𝑏₁₁+𝜎+𝜎𝜎²+𝑏₂₂𝜎−𝑏₂₄𝑏₄₂,(21) solving for the roots we get the following eigenvalues: 𝜎𝜎₁=−𝑏₁₁,(22)₂=−𝑏²²−𝑏₂₂²+4𝑏₄₂𝑏₂₄𝜎2,(23)₃=−𝑏²²+𝑏₂₂²+4𝑏₄₂𝑏₂₄2.(24)

The lower triangular −𝐁_𝟐 gives us our final two eigenvalues: 𝜎𝜎₄=−𝑏₃₃,(25)₅=−𝑏₅₅.(26)

ODE stability requires real parts of the 𝜎₁,…,𝜎₅ to be negative. In the next remark, stability criteria are formulated in terms of the parameters of the model.

Remark 4. Under the model assumptions we have −𝑏₁₁<0−𝑏₂₂<0,𝑏₄₂<0,−𝑏₅₅<0.(27) Therefore, 𝜎₁<0, 𝜎₂<0, and 𝜎₅<0. Next, if 𝑏₃₃=𝑎₂1−2𝐹𝑓^𝑒₀+𝑎¹²𝑐_𝑒−𝑎₃>0,(28) then 𝜎₄<0. Finally, because 𝑏₄₂<0, real part of 𝜎₃ is negative if and only if 𝑏₂₄=−𝑓₂𝜆₃−𝑓₃𝜆₂𝑐_𝑒>0.(29) Assumptions in (28) and (29) have clear biological interpretation.
Condition 𝑏₃₃>0 requires 𝑎²1−2𝐹𝑓^𝑒₀+𝑎¹²𝑐_𝑒<𝑎₃,(30) suggesting the need for the logistic growth of fibroblasts combined with the direct proliferation resulting from the presence of chemoattractants to be overcome by the death rate of fibroblasts.
Condition 𝑏₂₄>0 requires 𝑓₃𝜆₂𝑐_𝑒>𝑓₂𝜆₃,(31) suggesting that stability is aided when the percentage of regulatory macrophages out-weighs the percentage of inflammatory macrophages.
Note that from a mathematical point of view, conditions in the form of a strict inequalities imply a stronger property of the solution, namely asymptotic stability of the equilibrium. Lyapunov stability follows from the less restrictive condition with nonstrict inequalities.

4. ODE Linear Stability Implies PDE Linear Stability

Since the interior equilibrium solution represents the inflammatory state, one of the more biologically relevant questions is whether some modifications of conditions can cause the reactions to be away from the ill state and return to healthy state. Typically, the competition between diffusion and chemotaxis can aid the instability by creating spatial disturbance. One of the surprising findings for this system, however, is that if the equilibrium is stable by pure reactions, then it is stable for the whole reaction-diffusion-chemotactic system.

To start, we let 𝑣(𝑥,𝑡)=𝑒^𝜎𝑡𝜙_𝜇𝑛(𝑥)𝑢₁,…,𝑢₅(32) to be a vector with unknown five components and function 𝜙_𝑛(𝑥) to be the nth eigenfunction for Laplace equation with respect to Neumann boundary conditions: 𝜕𝜙_𝜇𝑛Δ𝜙_𝑛(𝑥)=−𝜇_𝑛𝜙_𝜇𝑛insidedomain,(33)𝜕𝑛=0ontheboundaryofthedomain.(34) Let us for simplicity assume that the domain is convex such that 𝜇_𝑛≥0 for any 𝑛∈ℕ is an eigenvalue for the eigenvalue problem, and 𝜙_𝜇𝑛 is its corresponding eigenfunction. We will drop the subscripts 𝑛 in the text below. Substituting the function 𝑣(𝑥,𝑡) into equation one can get 𝜎𝜎𝑢𝑢₁₂=−𝐷_𝑑𝜇𝑢₁−𝑏₁₁𝑢₁−𝑏₁₄𝑢₄,𝜎𝑢₃=−𝐷_𝑐𝜇𝑢₂−𝑏₂₁𝑢₁−𝑏₂₂𝑢₂−𝑏₂₄𝑢₄,𝜎𝑢₄=−𝐷_𝑓𝜇𝑢₃+𝜒_𝑓𝜇𝑢₂−𝑏₃₂𝑢₂−𝑏₃₃𝑢₃−𝑏₃₄𝑢₄,𝜎𝑢₅=−𝐷_𝑚𝜇𝑢₄+𝜒_𝑚𝜇𝑢₂−𝑏₄₂𝑢₂−𝑏₄₄𝑢₄,=−𝐷_𝑒𝜇𝑢₅+𝜒_𝑒1𝜇𝑢₃+𝜒_𝑒2𝜇𝑢₂−𝑏₅₃𝑢₃−𝑏₅₅𝑢₅,(35) or 𝑏𝜎+𝐷_𝑑𝜇+𝑏₁₁𝑢₁+𝑏₁₄𝑢₄=0,₂₁𝜎+𝐷_𝑓𝑢₁+𝜎+𝐷_𝑐𝜇+𝑏₂₂𝑢₂+𝑏₂₄𝑢₄=0,−𝜒𝜇+𝑏₃₃𝑢₃−𝜒_𝑓𝜇−𝑏₃₂𝑢₂+𝑏₃₄𝑢₄=0,_𝑚−𝜒_𝑒2𝜇−𝑏₄₂𝑢₂+𝜎+𝐷_𝑚𝜇+𝑏₄₄𝑢₄=0,𝜇𝑢₂−𝜒_𝑒1𝜇−𝑏₅₃𝑢₃+𝜎+𝐷_𝑒𝜇+𝑏₅₅𝑢₅=0.(36)

Then in matrix form it takes a form 𝐴(𝜎)𝑢=0,(37) with matrix 𝐴 defined as follows:⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝𝜎+𝐷𝑏_𝑑𝜇₂₁+𝑏₁₁00𝑏₁₄00−𝜎+𝜒𝐷_𝑐𝑓𝜇+𝑏₂₂0𝑏₂₄00−𝜒_𝑚𝜇−𝑏₃₂𝜎+𝐷_𝑓𝜇+𝑏₃₃𝑏₃₄00−𝜒𝜇_𝑒2−𝑏₄₂0𝜎+𝐷_𝑚𝜇+𝑏₄₄0𝜇−𝜒_𝑒1𝜇−𝑏₅₃0𝜎+𝐷_𝑒𝜇+𝑏₅₅⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠.(38)

Below, we will show that if the real part of all eigenvalues of matrix 𝐁 is negative (corresponding ODE system is stable), then nontrivial solutions of (37) with parameter 𝜎 having negative real part exist. It is not difficult to see that the determinant of the matrix 𝐴 has aform as follows: 𝑃(𝜎)=𝜎+𝐷_𝑒𝜇+𝑏₅₅𝜎+𝐷_𝑓𝜇+𝑏₃₃det𝐵₁.(39)

Here, 𝐵₁ is a matrix associated to debris 𝑢₁, chemotaxis 𝑢₂, and macrophages 𝑢₄ parameters only: ⎛⎜⎜⎜⎜⎝𝜎+𝐷_𝑑𝜇+𝑏₁₁0𝑏₁₄𝑏₂₁𝜎+𝐷_𝑐𝜇+𝑏₂₂𝑏₂₄0−𝜒_𝑚𝜇−𝑏₄₂𝜎+𝐷_𝑚𝜇+𝑏₄₄⎞⎟⎟⎟⎟⎠.(40)

Under the assumptions that the ODE part without diffusion is asymptotically stable, coefficients 𝑏₅₅ and 𝑏₃₃ should satisfy inequalities 𝑏₄₄=𝑏₁₄=0, 𝑏₅₅>0 and 𝑏₃₃<0.

We rearrange the matrix into a 𝑢₄,𝑢₁,𝑢₂ order so that it is similar to the one addressed previously in the ODE stability analysis. Now, det𝐁_𝟏+𝜎𝐈=⎡⎢⎢⎢⎢⎣𝜎+𝐷^𝑚0𝜇𝑏₁₁0𝑏₄₂−𝜒_𝑚𝜇𝑏+𝐷_𝑑𝜇+𝜎0₂₄=𝜎+𝐷_𝑚𝑏₂₁𝑏₂₂+𝐷_𝑐𝜇+𝜎⎤⎥⎥⎥⎥⎦−𝑏₂₄𝜇b₁₁+𝐷_𝑑𝜇+𝜎𝑏₂₂+𝐷_𝑐𝜇+𝜎=𝑏𝑏₄₂−𝜒_𝑚𝜇𝑏₁₁+𝐷_𝑑𝜇+𝜎₁₁+𝐷_𝑚+𝐷_𝑑𝜇+𝜎𝜎²+𝑏₂₂+𝐷_𝑐𝜇+𝐷_𝑚𝜇𝜎𝜇𝑏₂₂+𝐷_𝑐𝜇−𝑏₂₄𝑏₄₂−𝜒_𝑚𝜇,(41) solving for the roots we get the following eigenvalues: 𝜎𝜎₁=−𝑏₁₁−𝐷_𝑑𝜇,₂=−𝑏²²+𝐷_𝑐𝜇+𝐷_𝑚𝜇−𝑏₂₂+𝐷_𝑐𝜇+𝐷_𝑚𝜇²+4𝜖𝜎2,(42)₃=−𝑏²²+𝐷_𝑐𝜇+𝐷_𝑚𝜇+𝑏₂₂𝐷_𝑐𝜇+𝐷_𝑚𝜇²+4𝜖2,(43) here 𝜖=𝑏₄₂−𝜒_𝑚𝜇𝑏₂₄−𝐷_𝑚𝜇𝑏₂₂+𝐷_𝑐𝜇.(44) The other two eigenvalues are 𝜎𝜎₄₅=−𝑏₃₃−𝜇𝐷_𝑓,=−𝑏₅₅−𝜇𝐷_𝑒.(45)

In the forthcoming remark, explicit representations for all possible 𝜎’s are explored for direct comparison between conditions of the stability of the linearized PDE (8) and ODE (15) systems.

Remark 5. Similarly to criteria for ODE the stability for PDE, requires that real parts of the all 𝜎's to be negative. Under the natural constraints on the parameters of our original model 𝑏₁₁, 𝑏₅₅, 𝑏₂₂, and 𝑏₄₂ (see Remark 4) we already have 𝜎₁<0,𝜎₄<0 and 𝜎₅<0. Therefore, our criteria for PDE stability reduce to conditions as follows: 𝑏₂₂+𝐷_𝑐𝜇+𝐷_𝑚𝜇>0,𝜖<0.(46) It is obvious to see that if both inequalities hold, then 𝜎₂and𝜎₃ are negative. Since stability of the ODE system forces 𝑏₂₄>0 and 𝑏₄₂<0, these two inequalities for PDE stability hold for any 𝜒_𝑚>0, 𝐷_𝑚>0,𝐷_𝑐>0, 𝜇>0.
From the above arguments it follows that if the ODE system is stable, then 𝑣(𝑥,𝑡) are converging to zero as time goes to infinity for any eigenfunction 𝜙_𝑛. Therefore, since the 𝜙_𝑛(𝑥) is complete in 𝐿₂ space, one can conclude that the stability of the linearized PDE system (8) in 𝐿₂ space follows from the stability of the ODE system (15).
As expected, the ODE stability and PDE stability are different. Let 𝐷_𝑚=𝐷_𝑐𝜒_𝑚=0, then the first 5 eigenvalues of the PDE and ODE have the same sign. By definition of our original model 𝜎₁, 𝜎₂, and 𝜎₅ are all negative. Assume 𝑏₃₃>0 (in some sense reactive terms has stabilizing effect, with respect to 𝑈_𝑒), then 𝜎₄<0. However, now if one lets 𝑓₂𝜆₃>𝑓₃𝜆₂𝑐_𝑒, which means that inflammatory macrophages dominate the regulatory macrophages, then 𝑏₂₄<0 causing 𝜎₃>0, and consequently the ODE system (15) is unstable. For the same set of the coefficients 𝑏's and given 𝜇>0, it is not difficult to find sufficient condition on 𝐷_𝑚, 𝐷_𝑐, and 𝜒_𝑚 such that 𝜎₃<0, which guarantee stability of the equilibrium state 𝑈_𝑒. For example, any set with the same coefficients 𝑏’s with 𝐷_𝑚𝐷_𝑐>𝑏₄₂−𝜒_𝜇𝑏₂₄/𝜇(47) will have a real part of the 𝜎₃<0 and consequently the solution of the corresponding IBVP with initial function to be 𝜙_𝜇(𝑥)(𝑢₁,…,𝑢₅) will be vanishing at time infinity. Condition (47) contains the following pattern in the biological interpretation. Assume that inflammatory macrophages dominate the regulatory macrophages and are characterized by the coefficient 𝑏₂₄=−(𝑓₂𝜆₃−𝑓₃𝜆₂𝑐_𝑒)<0. Then for any given value 𝑏₂₄ if mobility of the macrophages and diffusion of the chemoattractant is high enough in comparison to the coefficient 𝑏₂₄, then 𝑈_𝑒 is stable for the class of perturbation which corresponds to eigenfunction 𝜙_𝜇. In less strict wording, the system can be cleaned from dead cells by high “mobility/diffusivity” of the macrophages with respect to chemoattractant. This indicates vital impact of the key parameters 𝐷_𝑚, 𝐷_𝑐, and 𝜒_𝑚 on “inflammatory” behavior both in space and in time of the system perturbed from equilibrium.
Obtained conclusion depends on 𝜇 and can be applied only if initial data is proportional to 𝜙_𝜇. If in the Fourier extension of the initial data all coefficients are nonzero, then the sufficient condition for stability is the same as for ODE system.
In the next section, we will analyze conditional stability of the IBVP for (8) under assumption that 𝑣(𝑥,𝑡₀) has zero average:∫𝑣(𝑥,𝑡₀)𝑑𝑥=0. We will derive conditions on the coefficient of the system (8) such that the 𝐿₂ norm of the solution is bounded by the 𝐿₂ norm of the initial data, or it converges to zero at time infinity depending on the conditions on coefficients. Those conditions will depend only on coefficients of the model and Poincare constant (𝐶_𝑝), which in turn depends only on the geometry of the domain. We will also show that there exists a specific initial distribution such that the corresponding IBVP solution is vanishing at time infinity while the corresponding solution of the ODE is unbounded at time infinity.

5. Stability of Equilibrium in the Linearized PDE System without ODE Stability

Let us rewrite the linearized system (8) as follows: 𝜕𝑑𝜕𝑡=𝐷^𝑑𝜆∇²𝑑−𝑓₀𝜆₁𝑚_𝑒𝑑−𝑏₁₄𝑚,(48)₁𝜕𝑐𝜕𝑡=𝜆¹𝜆𝐷_𝑐∇²𝑐−𝑏₂₁𝜆₁𝑑−𝑏₂₂𝜆₁𝑐−𝑏₂₄𝜆₁𝑚,(49)₁𝜕𝑓𝜕𝑡=𝐷^𝑓−𝑏𝜆₁₃₂∇²𝑓−𝜒_𝑓𝜆₁∇²𝑐𝜆𝜆₁𝑐−𝑏₃₃𝜆₁𝑓−𝑏₃₄𝜆₁𝑚,(50)₁𝜕𝑚𝜕𝑡=𝐷^𝑚−𝑏𝜆₁₄₂∇²𝑚−𝜒_𝑚𝜆₁∇²𝑐𝜆𝜆₁𝑐−𝑏₄₄𝜆₁𝑚,(51)₁𝜕𝑒𝜕𝑡=𝐷^𝑒−𝜒𝜆_1𝑒2∇²𝑒−𝜒_𝑒1𝜆₁∇²𝑓𝜆₁∇²𝑐−𝑏₅₃𝜆₁𝑓−𝑏₅₅𝜆₁𝑒.(52)

Next multiplying equations (48) by 𝑑, (49) by 𝑐, (50) by 𝑓, (51) by 𝑚, and (52) by 𝑒 correspondingly and integrating by parts, one can easily get 12𝜕𝜕𝑡𝑑²𝜆=−𝐷_𝑑(∇𝑑)²−𝑓₀𝜆₁𝑚_𝑒𝑑²−𝑏₁₄𝑚𝑑,(53)₁2𝜕𝜕𝑡𝑐²𝜆=−𝜆₁𝐷_𝑐(∇𝑐)²−𝑏₂₁𝜆₁𝑑𝑐−𝑏₂₂𝜆₁𝑐²−𝑏₂₄𝜆₁𝑚𝑐,(54)₁2𝜕𝜕𝑡𝑓²=−−𝑏₃₂𝐷_𝑓𝜆₁(∇𝑓)²+Φ(𝑐,𝑓)𝜆𝜆₁𝑐𝑓−𝑏₃₃𝜆₁𝑓²−𝑏₃₄𝜆₁𝑚𝑓,(55)₁2𝜕𝜕𝑡𝑚²=−−𝑏𝐷_𝑚42𝜆₁(∇𝑚)²+Φ(𝑐,𝑚)𝜆𝜆₁𝑐𝑚−𝑏₄₄𝜆₁𝑚²,(56)₁2𝜕𝜕𝑡𝑒²=+−Φ(𝐷𝑐,_𝑒𝜆𝑒)₁(−∇𝑏𝑒)²₅₃+Φ(𝑓,𝑒)𝜆₁𝑓𝑒−𝑏₅₅𝜆₁𝑒².(57) Here, Φ(𝑓,𝑒)=𝜒_𝑒1𝜆₁∇𝑓∇𝑒, Φ(𝑐,𝑓)=𝜒_𝑓𝜆₁∇𝑐∇𝑓, Φ(𝑐,𝑚)=𝜒_𝑚𝜆₁∇𝑐∇𝑚, Φ(𝑐,𝑒)=𝜒_𝑒2𝜆₁∇𝑐∇𝑒.

Adding LHS and RHS of the equations above: (53)+(54)+(55)+(56)+(57) and applying the Poincare inequality to the terms ∫(∇𝑢)²𝑑𝑥 such that for 𝐶_𝑝=𝐶_𝑝(Ω)>0, 𝐶_𝑝_Ω𝑢²𝑑𝑥≤_Ω(∇𝑢)²𝑑𝑥+_Ω𝑢𝑑𝑥²,(58) one can easily get 12𝑑²+𝜆₁𝑐²+𝑓²+𝑚²+𝑒²_𝑡≤−+−++𝐶𝐵𝐵[𝐵[((𝐵∇𝑐(,∇(𝑒𝑓𝑑,𝑐𝑑,)∇,𝑚∇+𝑐₂)𝑚)𝐵++)(+𝐵𝑓𝐵,𝐵((𝑚𝑐∇,()𝑑∇𝑒+,)𝑐∇+,𝐵∇𝑓𝐵(𝑓)𝑓](,)𝑐𝑒,)𝑚])+𝑐²+𝑓²+𝑚²+𝑒²,(59) where the bilinear forms are 𝐵(𝑐,𝑑)=𝜆₁𝐵(𝑑,𝑚)=0,(60)𝐵(𝑐,𝑚16)𝐷=^𝑐𝐶𝜆^𝑝₁+𝑏₂₂𝑐²+𝑏₂₁𝑑𝑐+𝑓₀+𝐷_𝑑𝐶_𝑝𝑑²,(61)+61𝐷1_𝑐𝐶_𝑝+𝑏₂₂𝑐²+𝑏_2,4𝑐𝑚3𝐷^𝑚𝐵(𝑓,𝑚)=𝜆₁𝐶_𝑝+𝑏₄₄𝑚²(62)𝐵(𝑐,𝑓)=𝜆41₁𝐷^𝑓𝐶_𝑝+𝑏₃₃𝑓²+𝑏₃₄𝑐𝑓+31𝐷^𝑚𝐶_𝑝𝑚²,(63)𝐵(∇𝑐,∇𝑚)=𝜆₁16𝐷𝐵^𝑐𝐶(𝑓^𝑝,𝑐𝑒²)+=𝑏0,³²𝑐𝑓+31𝐷^𝑓𝐶_𝑝𝑓²((6645))𝐵(∇𝑐,∇𝑓)=𝜆₁16𝐷^𝑐(∇𝑐)²−𝜒_𝑚∇𝑐∇𝑚+13𝐷^𝑚(∇𝑚)²,(66)𝐵(∇𝑓,∇𝑒)=𝜆₁61𝐷^𝑐(∇𝑐)²−𝜒_𝑓∇𝑐∇𝑓+14𝐷^𝑓(∇𝑓)²,(67)𝐵(∇𝑐,∇𝑒)=𝜆₁41𝐷^𝑓(∇𝑓)²−𝜒_𝑒2∇𝑓∇𝑒+12𝐷^𝑒(∇𝑒)²,(68)16𝐷^𝑐(∇𝑐)²−𝜒_𝑒1∇𝑐∇𝑒+12𝐷^𝑒(∇𝑒)².(69)

Imposing assumptions that all bilinear forms above are positively defined, one can then conclude that the system is stable. Below, we formulate a sufficient condition for the solution to be stable in 𝐿₂ space. The formulation of the assumptions is presented in terms of the parameters of the original system where biological meanings are more evident.

Condition 1. If 16𝐷^𝑐𝐶_𝑝+𝑓₃𝜆₂𝑚_𝑒+𝑓₄𝑓⁰+𝐷_𝑑𝐶_𝑝^1/2≥12𝑓¹,(70) then 𝐵(𝑐,𝑑)≥0.

Condition 2. If 14𝐷^𝑓𝐶_𝑝−𝑎₂1−2𝐹𝑓^𝑒₀+𝑎¹²𝑐_𝑒−𝑎₃13𝐷^𝑚𝐶_𝑝^1/2≥12𝑎¹𝜆₁,(71) then 𝐵(𝑓,𝑚)≥0.

Taking into account actual values for equilibriums 𝑐_𝑒 and 𝑓_𝑒 of the inflammatory equilibrium, one can reduce (71) to an inequality, which is easy to interpret.

Namely, assume that 41𝐷^𝑓𝐷_𝑚𝐶_𝑝+√𝐴𝐷_𝑚^1/2≥12𝑎¹𝜆₁,(72) then 𝐵(𝑓,𝑚)≥0. From the previously mentioned, 𝐴=𝑎₂−𝑎₃+𝑎₁₂𝑎𝑎⁰₁₁²+4𝑎²𝐹₀𝑎¹𝜆₁𝑓₄𝑓₀𝜆₁𝑎₀−𝑎₁₁𝑓₁𝑓₀𝜆₃𝑓₀𝜆₁𝑓₂𝜆₃𝑎₁₁−𝑓₃𝑎₀𝜆₂.(73) Due to the assumption (12), parameter 𝐴 is well defined for all values of the coefficients of the original model. Biological meaning of constraint (12) was explained in Remark 3, and it is necessary for the existence of the inflammatory equilibrium. What we want to point out here is that for any set of the parameters there exist large enough diffusive constants 𝐷_𝑚 and 𝐷_𝑓 that inequality (72) holds, and consequently bilinear form 𝐵(𝑓,𝑚)≥0.

Condition 3. If 16𝐷^𝑐𝐶_𝑝+𝑓₃𝜆₂𝑚_𝑒+𝑓₄13𝐷^𝑚𝐶_𝑝+𝑎₁₁𝑚_𝑒−𝑎₀²≥12𝑓²𝜆₃−𝑓₃𝜆₂𝑐_𝑒,(74) then 𝐵(𝑐,𝑚)≥0. For well posedness of the RHS in inequality (74), assume that 31=𝐷1_𝑚𝐶_𝑝+𝑎₁₁𝑚_𝑒−𝑎₀3𝐷^𝑚𝐶_𝑝+𝑎₁₁𝑓₄𝑓₀𝜆₁𝑎₀−𝑎₁₁𝑓₁𝑓₀𝜆₃𝑓₀𝜆₁𝑓₂𝜆₃𝑎₁₁−𝑓₃𝑎₀𝜆₂−𝑎⁰≥0.(75) We rewrite the above inequality in terms of the parameters of the original model to point out that for any given set of the parameters, there exists big enough coefficient 𝐷_𝑚, characterizing macrophages mobility, such that inequality (74) holds.

Condition 4. If 61𝐷^𝑐𝐶_𝑝14𝐷^𝑓𝐶_𝑝^1/2≥12𝑎¹²𝑓_𝑒,(76) then 𝐵(𝑐,𝑓)≥0.

Condition 5. If 61𝐷^𝑐13𝐷^𝑚^1/2≥12𝜒^𝑚,(77) then 𝐵(∇𝑐,∇𝑚)≥0.

Condition 6. If 16𝐷^𝑐14𝐷^𝑓^1/2≥12𝜒^𝑓,(78) then 𝐵(∇𝑐,∇𝑓)≥0.

Condition 7. If 41𝐷^𝑓12𝐷^𝑒^1/2≥12𝜒^𝑒2,(79) then 𝐷(∇𝑓,∇𝑒)≥0.

Condition 8. If 61𝐷^𝑐12𝐷^𝑒^1/2≥12𝜒^𝑒1,(80) then 𝐷(∇𝑐,∇𝑒)≥0.

We now assume that ∫ for all five components 𝑑(𝑥,0),𝑐(𝑥,0),𝑓(𝑥,0),𝑚(𝑥,0), and 𝑓(𝑥,0) is equal to 0 (initial data are orthogonal to 1. Then due to no-flux Neumann condition on the boundary for all times, _𝑈𝑑=_𝑈𝑐=_𝑈𝑓=_𝑈𝑚=_𝑈𝑒=0.(81) Therefore, the above Conditions (1–8) guarantee Lyapunov stability of the linearized system. If further for the same class of initial data we in addition assume strict inequalities in (70)–(76), then system will be asymptotically stable, and 𝐿₂ norm of the solution will exponentially converge to zero as time goes to infinity.

Here, we do not assume the ODE stability conditions of the equilibrium in this section. It will be easy to construct a specially inhomogeneous solution of the initial-boundary value problem (IBVP) so that the solution of corresponding ODE for 𝑉=∫𝑣(𝑥,𝑡)𝑑𝑥 is identically zero, where the PDE solutions can be either stable or unstable by adjusting certain parameters. Indeed, let the domain be a segment [0,𝜋] and as in (32), with 𝜙=cos𝑥. Then, in as Section 4, in order for 𝑣(𝑥,𝑡) to be a solution of corresponding IBVP it is necessary and sufficient that 𝜎 to be a root of the characteristic polynomial equation 𝑃(𝜎)=0 in (39). To see Conditions (1–8) are essential, we show an example of the system with: (1) Conditions (1–8) are all met, and (2)𝑃(𝜎) has a positive root in (39). For selected domain, assume Poincare constant 𝐶_𝑝=1. Assume that all coefficients are such that inequalities in all constraints except inequities in constraints Conditions 3 and 5 are satisfied. Let 𝑏₂₂≥4/5𝐷_𝑐, 𝑎₁₁𝑚_𝑒≥𝑎₀, and 0>𝑏₂₄≥−(𝐷_𝑐𝐷_𝑚/20)^1/2. Obviously for these set of the parameter Condition 3 satisfies. Then if √𝐷_𝑐𝐷_𝑚/60≥𝜒_𝑚 then Condition 5 holds and consequently 𝑣(𝑥,𝑡)→0 as 𝑡→∞. Furthermore, it is not difficult to see that if 𝑏₂₂=4/5𝐷_𝑐, and 𝑏₂₄=−(𝐷_𝑐𝐷_𝑚/20)^1/2 then in (42) is positive provided 𝜒_𝑚−𝑏₄₂𝐷^𝑐𝐷_𝑚20^1/2−95𝐷^𝑚𝐷_𝑐>0.(82) Inequality in (82) holds if 𝜒_𝑚>√102𝐷_𝑐𝐷_𝑚(83) Consequently, Condition (83) holds then ‖𝑣(𝑥,𝑡)‖_𝐿2→∞ as 𝑡→∞. Comparing stability in (77) and instability in (83), the conditions are optimal unto discrepancy in coefficients. In the next remark, we want to highlight the impact of the diffusive parameter and chemotactic coefficients on the stability of the inflammatory equilibrium 𝑈_𝑒.

Remark 6. In all above eight conditions inequalities hold for big enough values of diffusive coefficients 𝐷’s. This highlights the importance of the spatial distribution of the perturbation for the equilibrium. The major meaning of these condition is that for any set of the parameters if diffusivity coefficients are big enough then 𝑈_𝑒 is stable. Another key parameter, which characterizes the behavior of the spatial distribution of the system is the chemotactic coefficient. From the example above, one can see that if the chemotactic sensitivity coefficient 𝜒 is relatively bigger than the diffusivity characteristic of the process, then 𝑈_𝑒 is unstable. At the same time if it is relatively smaller, as in inequalities (77)–(80), then the inflammatory equilibrium is stable.

6. Conclusion and Discussion

To quantitatively study the processes governing inflammatory and fibrotic reactions against foreign bodies, we have built a mathematical model with the capability to predict the trends of macrophage migration, ECM production, and chemoattractant regulation by macrophages in these fibrotic reactions. The initiations of reactions are digestions of debris which are the natural responses of the immune system to damaged cells and tissues due to the implantation process. Our model is built based principally on biochemical mechanisms, and it has served its purpose in providing trends of reactions. The model is expressed by a system of partial differential equations with no flux boundary conditions.

We have considered an equilibrium state of the system and its stability conditions. We have provided a mathematical proof that when this equilibrium is stable in the corresponding ODEs, then it is also stable for the full system in 𝐿²(Ω). However, a system with a parameter set can be conditionally stable in the PDE sense when its ODE system is not necessarily stable. We provided some exclusive conditions for this to happen. These conditions correspond with feasible biological conditions, where the percentage of regulatory macrophages dominates that of the inflammatory macrophages.

We mention here that the system has infinitely many equilibria, all except for one containing at least one free parameter in it. The one under discussion here is called the interior equilibrium as it has 5 nonzero components. This particular equilibrium corresponds to an inflammatory state of the healing process, whose instability is an indicator of three possible dynamics: (1) best case scenario, returning to the healthy state; (2) uncertain development, transition to another “abnormal” equilibrium; (3) acute inflammatory response (worst case scenario), perturbations tend to infinity.

Our main mathematical result indicates that the inflammatory state’s stability mainly depends on the reaction dynamics and even that small spatial diffusion and big chemotaxis cannot destabilize the equilibrium which is stable in the reaction-only system. However, if the equilibrium is unstable by its reaction-only system, then spatial diffusion over chemotactic effects can help to stabilize the equilibrium if the initial perturbation is subjected to specific constraints. We did not discuss other equilibrium states due to the length of the paper, but there is no mathematical difficultly in accomplishing these tasks.

Acknowledgments

This work is supported by the National Institutes of Health Grant no. 1R01EB007271-01A2. The research of this paper was supported by the NSF grant DMS-0908177.