
==== Front
J Math Biol
J Math Biol
Journal of Mathematical Biology
0303-6812
1432-1416
Springer Berlin Heidelberg Berlin/Heidelberg

39012511
2126
10.1007/s00285-024-02126-4
Article
Building up a model family for inflammations
http://orcid.org/0000-0003-1442-1474
Reisch Cordula c.reisch@tu-braunschweig.de

1
Nickel Sandra 2
Tautenhahn Hans-Michael 2
1 https://ror.org/010nsgg66 grid.6738.a 0000 0001 1090 0254 Institute for Partial Differential Equations, Technische Universität Braunschweig, Universitätsplatz 2, 38106 Braunschweig, Germany
2 grid.411339.d 0000 0000 8517 9062 Clinic for Visceral, Transplantation, Thoracic and Vascular Surgery, Leipzig University Hospital, Liebigstrasse 20, 04103 Leipzig, Germany
16 7 2024
16 7 2024
2024
89 3 2920 12 2023
17 5 2024
7 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
The paper presents an approach for overcoming modeling problems of typical life science applications with partly unknown mechanisms and lacking quantitative data: A model family of reaction–diffusion equations is built up on a mesoscopic scale and uses classes of feasible functions for reaction and taxis terms. The classes are found by translating biological knowledge into mathematical conditions and the analysis of the models further constrains the classes. Numerical simulations allow comparing single models out of the model family with available qualitative information on the solutions from observations. The method provides insight into a hierarchical order of the mechanisms. The method is applied to the clinics for liver inflammation such as metabolic dysfunction-associated steatohepatitis or viral hepatitis where reasons for the chronification of disease are still unclear and time- and space-dependent data is unavailable.

Keywords

Mathematical modeling
Reaction–diffusion equations
Model hierarchy
Inflammation
Hepatitis
Mathematics Subject Classification

392-10
35K57
35Q92
http://dx.doi.org/10.13039/501100002347 Bundesministerium für Bildung und Forschung 031L0304C Tautenhahn Hans-Michael Deutsche Forschungsgemeinschaft465194077 428832822 Tautenhahn Hans-Michael Technische Universität Braunschweig (1042)Open Access funding enabled and organized by Projekt DEAL.

issue-copyright-statement© Springer-Verlag GmbH Germany, part of Springer Nature 2024
==== Body
pmcIntroduction

Modeling life science processes with only partly known mechanisms and actors is challenging. The challenge becomes even harder if there is only few data or only qualitative data available for validating the models. However, mathematical models for not fully understood complex problems provide insight into the relevance of mechanisms and allow testing hypotheses in-silico.

In this paper, a family of models for the dynamics of liver inflammation is formulated, analyzed, simulated, and considered in a clinical context.

Liver inflammations such as steatohepatitis or viral hepatitis often lead to chronic inflammation, which can subsequently result in serious diseases such as liver cirrhosis.

In case of viral hepatitis, the cause of the infection can now be eliminated very successfully using antiviral drugs so that the inflammation decreases accordingly. The trigger for steatohepatitis or autoimmune hepatitis is much more difficult to eliminate. Currently, immunomodulatory drugs are often used to combat sterile inflammation. These drugs possess a relevant side effect profile, but above all, the patient’s own immune system can respond far less well to pathogens. Mathematical models may help to bring order to the complex system of immune reactions and to better understand the course of the disease. In order to develop a suitable treatment strategy for such a complex immunological process as sterile inflammation, it is advantageous to know the exact mechanisms of the development and progression of inflammation. This would help to develop more specific drugs to control the inflammation.

A hierarchical order of the mechanisms is a priori unavailable because the mechanisms are coupled and can be isolated only in a few cases for single experiments. Models on a mesoscopic scale allow an abstraction from the unknown processes on the cell scale and allow comparing the modeling results to qualitative data like pathological images of the liver.

The model family describes the interactions and propagation of different influencing variables on inflammation. In particular, the following factors are currently considered in the model: macrophages (Kupffer cells), viruses, CD4+ T cells, cytotoxic T cells, and cytokines. Different approaches for known mechanisms are used and models of different complexity are analyzed and compared. As the focus lies on the model family, not every model will be discussed in all detail and not all mechanisms will be used in the models. The results for the model family provide insight into the hierarchical order of the mechanisms used.

Literature overview

The literature covers three topics: First, medical and biological observations on liver inflammations are given. Details on the known mechanism and involved cells follow in Sect. 2.2. Second, models for inflammations are presented, with a highlight on models for viral liver infections leading to inflammation. Third, some results on reaction–diffusion equations are given.

The leading causes of liver inflammation worldwide are viruses, metabolic diseases, autoimmune diseases, and alcohol. If left untreated for a long period of time, all of these factors can lead to cirrhosis of the liver, which is in fact life-threatening (Kazankov et al. 2019).

The course of viral liver infections begins with the infection itself, followed by a symptom-free incubation period. During an active phase, the immune system fights against the virus in the liver tissue. After the active phase, healing infection courses differ from chronic. While the amount of virus and the activity of the immune system reduces to zero in a healing infection course, the virus persists in the liver and the immune system remains active during a chronic infection course (Schiff et al. 2018). Inflammation is the process of an active immune system reacting to infected tissue. If the inflammation continues even for reduced viral loads, the inflammation is called chronic. Inflammation occurs in different tissues caused by virus, bacteria, or autoimmune diseases. This paper focuses on inflammation caused by viral liver infections like hepatitis B. The reasons for chronic infections are still partly unknown, compare (Bowen and Walker 2005; Thomas and Liang 2016), and data is rare. Qualitative data is given by pathological images of removed tissue (Kanel 2017). The pathological images show a partial inhomogeneous spread of T cells in the tissue parts considered. The amount of T cells is higher next to portal fields. During acute hepatitis B infections, up to 95% of the liver cells are infected with the virus (Kanel 2017; Schiff et al. 2018). During chronic phases, around one-third of the liver cells are infected (Kanel 2017; Schiff et al. 2018). These observations will be used as qualitative data in Sect. 3.3 for evaluating the modeling results.

Mathematical models for inflammations are available on different length scales and for different inflammation types. In (Aston 2018) models using ordinary differential equations are compared. Compartment models describe the total amount of free virus, healthy, and infected liver cells in the whole liver. Variations change the reaction functions or take delay into account, compare (Nangue et al. 2022). The same components are considered in Rezounenko (2018) but the diffusion of all three cell types are added. In this approach, even the liver cells are diffusing, which is not realistic. A more realistic model that is adapted to hepatitis B infections is presented in Tadmon et al. (2021). The modeled cells are three types of liver cells, free virus and B cells, and only the free virus is diffusing. This leads to a coupled system of ordinary and partial differential equations for which stationary states are analyzed and chronic infections are investigated.

The presented models describe the dynamics of the liver cells and the free virus. As available data is given by pathological images showing the T cells, the results are only partly comparable with the space-dependent images.

Mathematical models need to be spatially resolved for comparison with pathological images. Additionally, the T cells need to be modeled as these are the only cell types next to liver cells displayed in the pathological images. Reaction–diffusion equations are one possibility for modeling interacting cells in space and time. They allow a mathematical analysis and in some cases a prediction of the solution behavior.

Some results on reaction–diffusion equations are summarized here. Reaction–diffusion equations are used for modeling many applications like morphogenesis (Turing 1952), the spread of populations (Murray 2002) or chemical reactions (Schnakenberg 1979). Depending on the reaction functions and the diffusion parameters, the models show different types of solutions like traveling waves, (Turing) patterns, blow-ups, or relaxation (Perthame 2015). In some cases, these solutions can be proven and predicted analytically. In the light of chronic and healing infection courses of liver inflammations, two types of solutions are particularly interesting: Solutions tending towards zero describe healing infection courses and chronic infection courses are mathematically described by solutions with a tendency to stationary and spatially inhomogeneous distributions. As there are some asymptotic results on solutions tending towards spatially homogeneous states, c.f. Smoller (1994), the contradiction of these results can be used for gaining requirements for solutions not tending towards a homogeneous steady state, compare (Reisch and Langemann 2020).

Overview of this paper

The model family is formulated in Sect. 2. Biological information is translated into classes of feasible reaction functions (Sect. 2.2.1) and of taxis terms (Sect. 2.2.2). Sections 2.3 and 2.4 give examples of these classes. The model family is analyzed in Sect. 3 and requirements on the classes are formulated in Sect. 3.2.4. Section 3.3 presents numerical simulations and compares them with qualitative observations. Section 3.4 highlights the clinical impact of mathematical models in the context of inflammations, and a conclusion follows in Sect. 4.

Model family

A model family consists of models describing an application using the same modeling approach, but differing in the number or complexity of mechanisms included. First, we discuss the model type of reaction–diffusion equation and boundary conditions. Afterward, biological information on the dynamics of liver infections is given and is translated into conditions on function classes.

Reaction diffusion equations

The models have the general form1 q,t=F(t,x,q)+∇·D(x,q,∇q),

where t>0, x∈Ω⊂Rd and Ω is a model for the part of the liver considered. In what follows, the dimension of space is d∈{1,2,3}, depending on the modeling of the domain.

The vector-valued function q:[0,∞)×Ω→Rn gives the time and space-dependent amount of the acting cells and virus. The function F describes the reactions between the cells and the virus, as well as the growth or decay of any substances described. The term D(x,q,∇q) models any diffusive or chemotactic effects caused by gradients of substances or external influence.

The reaction–diffusion in Eq. (1) is completed to an initial value problem by q(0,x)=q0(x) for x∈Ω with q0(x)≥0 for all x∈Ω. Additionally, we assume zero-flux boundary conditions2 ∂q∂n=0forx∈∂Ω.

The homogeneous Neumann boundary conditions are an approximation: The liver consists of eight different segments which are separated by infeasible tissue, Schiff et al. (2018). There, zero flux conditions are reasonable because the exchange with the exterior of one segment happens via the blood vessels.

If Ω is a part of such a liver segment, homogeneous Neumann boundary conditions are less intuitive. Assuming that the liver tissue in Ω is homogeneous, periodic boundary conditions are a reasonable choice. Biologically, liver inflammations are stronger in some parts of the liver than in others. Therefore, the assumption of similar structures next to each other is not perfect.

As an approximation for nearly constant amounts of acting cells close to the boundary ∂Ω, zero flux conditions are acceptable and reflect the higher importance of local reactions compared to diffusion over the boundary.

Biological motivation

Viral hepatitis such as hepatitis B is named after the virus type causing an infection and leading to inflammation of the liver tissue. Inflammation is a consequence of the reaction of the immune system to the infection.

Hepatitis B infections show two typical infection courses. In healing infection courses, the virus and the inflammation vanish after an acute inflammation with a strong active phase of the immune system. Chronic infection courses have an active phase followed by a chronic phase with mild symptoms. In the chronic phase, virus and inflammation persist in the liver tissue. Chronic liver infections often lead to severe secondary diseases like liver cancer or cirrhosis.

The following is a brief summary of the cellular immune response. After infection, the virus spreads in the body and settles in the liver. There it multiplies and spreads throughout the liver. In this phase, the adaptive immune response begins to take effect. Dendritic cells recognize and process virus particles and present them to T cells. The T cells specify and develop to CD4+ or CD8+ T cells. Some CD4+ T cells become T helper cells, specializing in different types of T helper cells. One type of T helper cell enhances the cytokine concentration, which are proteins caused by the interaction of tissue cells and virus. Cytokines are important for cell signalling and lead the way for T cells towards the virus.

The CD8+ T cells named cytotoxic T cells trigger the cell death of infected tissue cells. Therewith, the tissue cell and the contained virus are destroyed. Inflammation is essential to neutralize the damaging stimulus and create the conditions for repair processes. The inflammation is caused only indirectly by the virus.

Next, the description of infection is split up into single mechanisms. From a model-theoretic perspective, this step is again a modeling process. We start with naming the actors of the dynamics. Let q1=q1(t,x) be the space and time-dependent amount of the virus. According, q2=q2(t,x) is the amount of T cells, which can be divided into T helper cells Th=Th(t,x) and cytotoxic T cells Tc=Tc(t,x). Cytokines q3=q3(t,x) are the third quantity.

The effects on each other are visualized by arrows or dependencies (∼). virus replicate and spread out: q1→q1↗

immune system reacts to virus by producing T cells: q1→q2↗

T helper cells produce and enhance cytokines: Th→q3↗

cytokines guide cytotoxic T cells: q3∼Tc

cytotoxic T cells kill virus: Tc→q1↘

in absence of virus, T cells reduce (as well natural decay): ¬q1→q2↘

random spreading of T helper cells

Fig. 1 Mechanisms during an inflammation. The variables are virus q1, cytokines q3 and T cells q2, divided into T helper cells Th and cytotoxic T cells Tc. Positive effects are displayed as +, negative effects as -, and movement effects as m. The spreading of the virus and T cells is not visualized. The dashed arrow shows the enhancement of mechanism M3 depending on q1

Figure 1 visualizes the mechanisms and their effects. There is a large loop effect: Mechanism M2 enhances the production of both T cell types according to the virus amount q1. The T helper cells Th produce, again depending on q1, cytokines. The cytokines q3 control the cytotoxic T cells Tc, which reduce the virus q1.

Reactions

The next step in building up a model family is the translation of the biological observations into mathematics. Figure 1 is an intermediate step for this process. There, the mechanisms are converted from a verbal description to a logical expression. Next, the logical expressions are translated into mathematical formulas, regarding more details of the verbal expression than in Fig. 1.

Therefore, we specify the reaction–diffusion system in (1). F1 is the first entry of the vector F and belongs to the partial differential equation for the virus q1. Analogously, F2 and F3 belong to q2 and q3. If the T cell population in the model is split up into the two subtypes Th and Tc, there are reaction functions FTh and FTc.

Every reaction function Fi consists of different single mechanisms fiMj splitting up the biological mechanisms occurring during an inflammation. The reaction function for the change of the virus is given by3 F1(t,x,q)=f1M1(t,x,q)+f1M5(t,x,q),

where the first reaction function f1M1 describes the growth and f1M5 describes the decay of the virus. The spreading of the virus, which is also part of mechanism M1, is described by a taxis term, which will be discussed in Sect. 2.2.2.

If the T cells are only considered as general T cells fulfilling all mentioned functions of T helper cells and cytotoxic T cells, then the biological mechanisms affecting q2 areF2(t,x,q)=f2M2(t,x,q)+f2M4(t,x,q)+f2M6(t,x,q).

If the T cell subtypes are split up into T helper cells Th and cytotoxic T cells Tc, then the mechanisms affect directly the subtypes instead of the general T cells, soFTh(t,x,q)=fThM2(t,x,q)+fThM6(t,x,q),FTc(t,x,q)=fTcM2(t,x,q)+fTcM4(t,x,q)+fTcM6(t,x,q).

So far, the reaction function for the cytokines consists of one mechanism,F3(t,x,q)=f3M3(t,x,q)

describing the production of cytokines by T helper cells in the presence of the virus.

The mechanisms translate into conditions on the partial derivatives of the functions fiMj. By considering first the translation into partial derivatives and not directly into certain functions, the classes of reaction functions include more variety. The deductive modeling process from biological mechanisms to a mathematical model family is highlighted by these finer steps.

In (Reisch and Langemann 2019a) the conditions∂F1∂q2<0,∂F2∂q1>0,∂F2∂q2<0

are presented, which are now investigated in more detail. We start with the mechanisms for the virus. The growth of the virus (M1) can be modeled by a function fulfilling4 ∂f1M1(t,x,q)∂q1≥0.

The equality with zero can be interpreted as a stopped growth of the virus, for example, due to an eliminated virus or due to saturation. The second mechanism affecting the virus is a reduction depending on the cytotoxic T cells. This mechanism can be expressed by the conditions on the partial derivatives5 ∂f1M5(t,x,q)∂Tc≤0or∂f1M5(t,x,q)∂q2≤0,

depending on the chosen variables.

Analogously, we find for the T cells6 ∂f2M2(t,x,q)∂q1≥0,∂f2M4(t,x,q)∂q3≥0,∂f2M6(t,x,q)∂q2≤0.

The production of T cells in mechanism M2 models the impact of the immune system on the virus. The more virus there is, the more T cells are produced. The mechanism M4 shows that the amount of T cells increases depending on the cytokines. In mechanism M6, the decay of T cells is modeled.

If the model considers a splitting up of the general T cells into T helper cells and cytotoxic T cells, then there are conditions on the mechanism functions for the special T cells. For the T helper cells, the mechanisms M2 and M6 are relevant and the conditions in (6) hold for fThM2 and fThM6. The mechanisms M2, M4 and M6 are relevant for the cytotoxic T cells and (6) holds for fTcM2, fThM4 and fTcM6.

The conditions for the reaction functions of the cytokines are7 ∂f3M3(t,x,q)∂Th≥0and∂f3M3(t,x,q)∂q1≥0,

where the partial derivative with respect to the T helper cells is replaced by the partial derivative with respect to the T cells q2 in the case of not splitting up into subtypes.

For all quantities, a natural decay will be discussed in Sect. 2.3. The natural decay is a result of different biological processes and a small amount of transport. For the virus, the natural decay can be interpreted as a reduced growth factor. Additionally, some processes have a natural saturation, for example, due to restricted resources.

The general conditions on the partial derivatives of the mechanism functions are one step of the modeling process, providing the substance for building a model family. In the next section, conditions on the taxis terms are presented. Afterward, particular reaction functions and taxis terms fulfilling the conditions are given and the effects of natural decay and saturation are discussed in more detail.

Taxis

The taxis are not fully covered by the description in Fig. 1 and are therefore explained in more detail here. The domain Ω is a part of the liver, not containing, for example, the lymph. Consequently, the details of the transport of T cells into the liver is not part of the model. The domain Ω contains a subdomain Θ where the inflow of T cells takes place. If Ω is a two-dimensional domain, Θ can be considered as a cut through a portal field of the liver structure.

The boundary conditions prevent an exchange of virus or T cells with the exterior. In this section, the taxis of any cells in the interior of Ω are described.

The virus q1 spreads out in the liver. Two different mechanisms cause the spreading (Goyal and Murray 2016). Cell-to-cell transmission leads to a spreading with a small medium step size, and the virus in one infected liver cell passes on to a neighbor cell after replication. The second mechanism is the diffusion of the virus through the extracellular space. The medium step size is typically larger than for the cell-to-cell transmission. The diffusion function in (1) for the q1-component therefore consists of two termsD1(x,q,∇q)=D1ctc(x,q,∇q)+D1ecs(x,q,∇q).

The cell-to-cell transmission can be interpreted as a diffusion process with a small diffusion parameter. This mechanism has the form8 D1ctc(x,q,∇q)=D1ctc(∇q1).

The spreading through the extracellular space might depend on the blood stream through the extracellular space and therefore on the direction and velocity of the blood. The results term would include this extracellular spread additionally to the virus amount, for example9 D1ecs(x,q,∇q)=D1ecs(x,q1,∇q1).

The cytotoxic T cells Tc follow a cytokine gradientDTcchem(x,q,∇q)=DTcchem(Tc,∇q3)

and the taxis depend on the cytotoxic T cells and the cytokine gradient.

The T helper cells spread by diffusion, so DThdiff(x,q,∇q)=DThdiff(Th).

For the general T cells without a division into subtypes, both effects combine and the dependencies on Th and Tc are replaced by a dependency on q2.

Cytokines are proteins important for cell signaling without self-induced movement. Due to natural processes, a small diffusion is possible. Additionally, the transport of the cytokines by the bloodstream through the extracellular tissue might happen. As a combination of those two mechanisms, the term for the cytokine q3 isD3(x,q,∇q)=D3diff(∇q3)+D3ecs(x,q3,∇q3).

We give examples for taxis terms in Sect. 2.4.

Examples of feasible reaction functions

The model in Eq. (1) is specialized for particular models of the model family. The functions are motivated by biology, see Sect. 2.2, and still variable enough to cover various scenarios.

Mechanism M1 describes the virus replication, and functions fulfilling Eq. (4) are10 f1M1(q1)=a1q1(unbounded),f1M1(q1)=a1q1(C1-q1)(bounded byC1),f1M1(q1)=a1q1(C1-q1)q1-εq1+κ(bounded byC1and Allee effect).

The three examples are extensions of each other. The first function describes an unbounded exponential growth of the virus. The second function is bounded from above by a constant capacity C1. The resulting growth is a logistic growth, starting exponentially and reaching saturation, which describes that all liver cells are infected.

The third example includes a strong Allee effect (Allee 1949), modeling a small vanishing virus amount by a negative function for q1<ε. This effect can be interpreted as a small immune reaction acting on tiny viral loads. A discussion of the differences in the dynamical system’s behavior due to the Allee effect can be found in Reisch and Langemann (2019a). The function with the Allee effect fulfills the condition in Eq. (4) only for values q1≥ε.

The mechanism M2 describes the production of T cells, both T helper cells and cytotoxic T cells if the virus is present. Possible mechanisms are11 f2M2(q1)=a2q1(local, unbounded),f2M2(q1,q2)=a2q1(C2-q2)(local, bounded byC2),f2M2(x,q1)=a2χΘ(x)∫Ωq1dx(global),f2M2(x,q1)=a2χΘ(x)(C2-q2)∫Ωq1dx(global, with saturation).

The first example depends solely on the local virus amount. T cells appear where the virus is. This effect is extended by an upper bound C2 for the T cells in the second example. The local mechanism neglects the finer liver structure with portal fields. T cells are produced outside of Ω and are transported via blood vessels into the liver. In the small-scale liver structure, the influx of T cells takes place in portal fields Θ⊂Ω.

The subdomains Θ are part of the third and fourth examples for mechanism functions. The functions are explicitly space dependent via a function χΘ(x) fulfilling12 ∫ΩχΘ(x)=1andχΘ(x)≥0forx∈Θ,=0forx∈Ω\Θ,

compare (Reisch and Langemann 2019a). The function χ can be a characteristic function for the subdomain Θ, or, if required, a smooth function with positive values and a local support in Θ.

The integral ∫Ωq1dx models the total viral load in Ω. The production of T cells is proportional to the total viral load in examples three and four in Eq. (11).

If the T cells are distinguished into T helper cells Th and cytotoxic T cells Tc, the proposed mechanisms are adapted by replacing q2 by one of the subtypes with individual capacities CTh, CTc. The requirement in Eq. (6) is fulfilled for all examples.

Mechanism M3 describes the production of cytokines by the virus and the effect of T helper cells increasing the cytokines. Examples of reaction functions are13 f3M3(q1)=a3q1(only dependent on virus),f3M3(Th)=a3Th(only dependent on T helper cells),f3M3(q1,Th)=a3q1Th(unbounded),f3M3(Th,q3)=a3Th(C3-q3)(dependent onThcells, bounded),f3M3(q1,Th,q3)=a3Thq1(C3-q3)(bounded),

where Th can be replaced by q2 if there is no division into subtypes of T cells. The examples differ in their dependencies on the involved cells and in the boundedness of cytokines at a certain place. Dependency only on the virus neglects the mechanism that T helper cells increase the cytokines, while functions independent of the virus neglect the mechanism that cytokines are produced by infected liver cells.

The conditions in Eq. (7) are fulfilled. The last two examples in Eq. (7) depend on q3 itself and give an upper bound. By only regarding mechanism M3, the constant C3 gives an upper bound. Further discussion on the boundedness of solutions follows in Sect. 3.2.

The examples for mechanism M3 in Eq. (13) can be ranked hierarchically. The more mechanisms are included, the higher the ranking is. The bottom level contains only the first two examples. The last example contains the most mechanisms and is therefore at the top of the ranking. Containing the most mechanisms is not always a feature of quality, as the mechanisms and their importance are unknown.

M4 describes the control of the movement of the cytotoxic T cells, c.f. Sect. 2.4.

The mechanism M5 describes the decay of the virus in dependence on cytotoxic T cells. Possible functions for this mechanism are14 f1M5(Tc)=-a5Tc(only dependent on cytotoxic T cells),f1M5(q1,Tc)=-a5q1Tc(dependent on virus and cytotoxic T cells).

The condition in Eq. (5) is fulfilled by both functions in (14).

Again, the dependency on the cytotoxic T cells can be replaced by a dependency on the general T cells. The first example in Eq. (14) does not preserve the non-negativity of q1. The analysis of the reaction functions follows in Sect. 3.2 in more detail.

The mechanism M6 describes that the T cell amount does not decrease if the virus is still present. It includes different effects: First, there is a decay of T cells. Second, this decay is smaller if there is more virus. The technical description of ’a lot of virus’ depends on the boundedness of the virus q1 by C1. Examples are15 f2M6(q2)=-a6q2(natural decay),f2M6(q1,q2)=-a6q2(C1-q1)(decay in absence of virus),

where the first example is included in the second becausef2M6(q1,q2)=-a6q2(C1-q1)=-a6C1q2+a6q1q2=C1f2M6+a6q1q2

yields. Both examples in Eq. (15) fulfill the condition in Eq. (6) on the partial derivative with respect to the T cells. The second example in Eq. (15) consists of a natural decay and an increase depending on the virus and the T cells. The increase can be interpreted as a benefit of the T cells, like in a predator–prey system or a host-pathogen system. Further details on this comparison are discussed in Reisch and Langemann (2019b). If the T cells are divided into T helper cells and cytotoxic T cells, mechanism M6 works in the same way by replacing q2 by Th or Tc, respectively.

Examples of the reaction functions of the different mechanisms are given in Eqs. (10), (11), (13), (14) and (15). The taxis describing mechanism M4 will be discussed in Sect. 2.4. Additionally, there is a natural decay for all quantities.

The virus is decaying in the presence of (cytotoxic) T cells, see mechanism M5 in Eq. (14). Additionally, there is a small natural decay, which is already included in the growth, covered by mechanism M1 in Eq. (10). The growth incorporating an Allee effect exhibits a negative value for a tiny amount of virus. This negative growth can be interpreted as a natural decay of a tiny amount of virus. Besides, a natural decay can be expressed by a smaller parameter a1. Consequently, the natural decay for the virus is already included in the example for mechanism M1 with the Allee effect. The other examples in Eq. (10) do not include the natural decay, interpreted as a local immune reaction, for small amounts of virus in the growth function. This results in different analytical properties, see Sect. 3.2.

The reaction function of mechanism M6 for the T cells contains a natural decay, too. In the case of a division of the T cells into the subtypes of T helper cells and cytotoxic T cells, the mechanism M6 is used in the same way for both subtypes.

So far, a natural decay for q3 like16 f3nd(q3)=-andq3,

is not contained in the mechanisms. A possible constant decay is17 f3nd=-andifq3>0,0ifq3=0,

resulting in an unsteady reaction function.

This section presents examples of reaction functions fulfilling the requirements from Sect. 2.2.1 and describing mathematically the biologically motivated mechanisms. Next, we give examples of the movements of the different quantities.

Examples of feasible taxis terms

The general reaction–diffusion equation in Eq. (1) has a taxis term of the form ∇·D(t,x,q,∇q). This term was divided into terms for every component in Sect. 2.2.2. In this section, the taxis terms for every component are specified and examples are given.

The virus q1 spreads out using cell-to-cell transmission and diffusion through the extracellular space. Possible dependencies for the two spreading mechanisms are given in Eqs. (8) and (9). Cell-to-cell transmission is a pure diffusive process, soD1ctc(∇q1)=d1ctc∇q1.

The spreading of the virus through the extracellular space can be modeled asD1ecs(∇q1)=d1ecs∇q1(homogeneous diffusion),D1ecs(x,∇q1)=d1ecsA(x)∇q1(space dependent diffusion),

where A(x)∈Rd×d is a positive definite, space-dependent matrix modeling the flux directions in the extracellular space at every point x∈Ω.

In more detailed modeling approaches, the mechanisms could be additionally time-dependent and display structural deformation on cell-scale caused by inflammation.

The taxis of T cells is type-dependent. Cytotoxic T cells are attracted by cytokines, whereas the taxis of T helper cells is dominated by diffusion. We start with the diffusion of the T helper cells, which is analogous to the diffusion of the virusDThdiff(∇Th)=dThdiff∇Th.

The motion of the cytotoxic T cells is dominated by chemotactic effects, which means that the cytotoxic T cells follow a cytokine gradient. Examples are18 DTcchem(∇q3)=dTcchem∇q3,DTcchem(Tc,∇q3)=dTcchemTc∇q3,

where the first example in Eq. (18) does not depend on the number of cytotoxic T cells and might lead therefore to a negative value Tc.

If the T cells are not divided up into the two subtypes, there might evolve conflicts with the dynamic of T cells and cytokines. These problems will be discussed later. An addition of the motion of cytotoxic T cells and T helper cells leads to a taxis termD2(q2,∇q2,∇q3)=d2diff∇q2+d2chemq2∇q3.

The cytokines are spreading by some diffusion D3(∇q3)=d3diff∇q3.

Next, different models are built by using the variety in the model family.

Analysis of the model family

The reaction functions and taxis terms in Sect. 2.3, and 2.4 are combined to models of the model family. Afterward, analytical results for the models are discussed. This analysis gives further requirements on the variety of possible reaction and taxis functions. Numerical simulations of acceptable models follow in Sect. 3.3.

Examples for models

The model family includes models for different sets of cell types. The largest set is the combination (q1,Th,Tc,q3). Another possible set of participating quantities is (q1,q2,q3). The smallest model with interactions consists of virus and T cells, (q1,q2) and neglects the cytokines. For certain analysis steps, different combinations or one-component models are meaningful as well, see Reisch and Langemann (2019a, 2020).

Remark 1

(Model 1) A large model in our model family for the four components virus q1, T helper cells Th, cytotoxic T cells Tc, and cytokines q3 is19 q1,t=a1q1(C1-q1)q1-εq1+κ-a5q1Tc+∇·[d1ctc∇q1+d1ecsA(x)∇q1],

20 Th,t=a2,hχΘ(x)(CTh-Th)∫Ωq1dx-a6,hTh(C1-q1)+∇·[dThdiff∇Th],

21 Tc,t=a2,cχΘ(x)(CTc-Tc)∫Ωq1dx-a6,cTc(C1-q1)-∇·[dTcchemTc∇q3],

22 q3,t=a3Thq1-andq3+∇·d3diff∇q3,

for x∈Ω and t>0. Let q=(q1,Th,Tc,q3)T. Boundary conditions are like in Eq. (2) zero flux conditions. Non-negative initial values complete the problem. All parameters have positive values. The matrix A is positive definite and χΘ fulfills Eq. (12).

Model 1 in Remark 1 uses all mechanisms and for every mechanism the example with the highest complexity in comparison to the others. Of course, there are many more detailed functions and models with many more components thinkable.

Next, a model with three interacting cell types is presented, compare (Reisch and Langemann 2019a).

Remark 2

(Model 2) A second model for inflammation contains the three components virus q1, T cells q2 and cytokines q3. The dynamics are given by23 q1,t=a1q1(C1-q1)q1-εq1+κ-a5q1q2+∇·[d1ctc∇q1],

24 q2,t=a2χΘ(x)∫Ωq1dx-a6q2(C1-q1)+∇·[d2diff∇q2-d2chemq2∇q3],

25 q3,t=a3q1-andq3+∇·[d3diff∇q3],

again completed with zero flux boundary conditions Eq. (2) and non-negative initial values. The parameters are positive and χΘ fulfills Eq. (12).

As a further reduction, the cytokines are not considered. The right circle in Fig. 1 containing the mechanisms of the cytokines is neglected. The resulting model only contains the dynamics of virus and T cells, compare (Kerl et al. 2012; Reisch and Langemann 2019a).

Remark 3

(Model 3) A two-component model of virus q1 and T cells q2 is26 q1,t=a1q1(C1-q1)q1-εq1+κ-a5q1q2+∇·[d1ctc∇q1],

27 q2,t=a2χΘ(x)∫Ωq1dx-a6q2(C1-q1)+∇·[d2diff∇q2].

Again, the parameters are positive, χΘ fulfills Eq. (12), the initial conditions are non-negative and there are zero flux boundary conditions.

Model 3 in Remark 3 uses different examples of reaction functions that have fewer effects than model 1 in Remark 1. This is not only a consequence of the reduction of modeled cell types but also depends on the reaction functions used.

It is not possible to decide a priori, that means before an analytical analysis and simulations, for the best model. This is caused by the uncertainty in choosing mechanisms, reaction and taxis functions. Further, it is a consequence of a missing definition for a ‘best’ model.

With these three different models in mind, we start the analysis of the models with as general reaction functions as possible.

Analytical results

The analysis of the reaction–diffusion equations covers some basic properties.

Non-negativity of solutions

The first property which will be discussed is the non-negativity of the solutions. As the solutions are interpreted as an amount of a certain cell type, negative values are meaningless in the light of application.

The reaction–diffusion equations consist of some reaction functions and divergence terms. If the divergence term only includes diffusive effects, the quantity influenced by the effect will not become negative. In those cases, it is sufficient to regard the reaction terms and check whether it allows a negative quantity. We discuss first the reaction functions and afterward the taxis.

The reaction function for the virus q1 consists of two mechanism functions for M1 and M5, see Eq. (3). The mechanism M1 depends on the virus, so f1M1=f1M1(t,x,q1). The explicit time and space dependency is only theoretical and not discussed.

The mechanism M5 depends on the amount of (cytotoxic) T cells, Tc or q2, and the amount of virus. The dependency reads f1M5=f1M5(t,x,q1,Tc), where Tc can be replaced by q2. Together, the reaction function is called F1(t,x,q1,Tc), and it specifies Eq. (3). A necessary condition for the non-negativity of q1 is28 F1(t,x,0,Tc)≥0.

A function modeling the growth of a virus should fulfill29 f1M1(t,x,0)=0.

All examples in Eq. (10) fulfill Eq. (29), so there is no increase in virus if no virus is present.

The mechanism M5 depends on the virus itself and on the amount of (cytotoxic) T cells. Additionally, the requirement Eq. (5) is a non-positive partial derivative with respect to the T cells. Therefore, f1M5 has negative values, and should be zero if q1=0. Otherwise, there would be an automated cell death without any virus. Again, the requirement in Eq. (28) transfers to a condition on f1M5(t,x,q1,q2) as30 f1M5(t,x,0,q2)=0.

Regarding the examples of reaction functions for mechanism M5 in Eq. (14), only the second example fulfills these requirements. The first example depends solely on the T cells and can lead to a negative amount of virus if q1=0.

Next, we regard the reaction functions of the T cells q2. The discussion holds as well for the subtypes T helper cells Th and cytotoxic T cells Tc. The reaction function consists of two mechanisms, M2 for the production and M6 for the decrease. The production depends on the virus and, if a saturation effect is considered, on the amount of T cells. So, f2M2=f2M2(t,x,q1,q2). The decay of T cells depends as well on the virus and on the amount of T cells. While the dependency on q1 is a fundamental requirement in mechanism M2, the dependency on q1 in M6 is optional, compare Eqs. (11) and (15). In total, the reaction function for the T cells and the subtypes is a function F2(t,x,q1,q2). As the function of M2 is non-negative, the requirement31 F2(t,x,q1,0)≥0

transfers directly to a requirement on the reaction function of mechanism M6, so32 f2M6(t,x,q1,0)≥0.

All examples in Eqs. (11), (15) fulfill Eq. (32) and are therefore, concerning the non-negativity of solutions, suitable functions. Of course, other functions are thinkable as long as the non-negativity requirements are fulfilled.

The reaction function F3 of the cytokines q3 consists of functions for the mechanism M3 and on the natural decay, compare Eqs. (13), (16), (17). Mechanism M3 describes the increase in the quantity of cytokines and can depend on the virus, the T (helper) cells, and the quantity of cytokines, so f3M3=f3M3(t,x,q1,q2,q3). The natural decay depends only on the quantity of cytokines, either directly like in Eq. (16), or indirectly like in Eq. (17). In total, the reaction function F3 of the cytokines depends on all cells, so the non-negativity requirement reads33 F3(t,x,q1,q2,0)≥0.

The function f3M3 describes an increase and therefore is non-negative. A function modeling the natural decay of cytokines should fulfill the requirement34 f3nd(t,x,0)=0

to model the natural decay. Both examples in Eqs. (16), and (17) fulfill the condition.

As seen, most of the conditions on the reaction functions for preserving non-negativity transfer directly into conditions on single reaction functions. Most of the introduced examples for reaction functions fulfill these conditions.

As a next step, the taxis terms are considered. We start with the diffusive movement, described by a term ∇·[didiff∇qi]=didiffΔqi.

The diffusive spreading is well known from the heat equation. It cannot result in a negative amount of qi if qi is non-negative for all x∈Ω. This is a consequence of the relaxation nature of diffusion and can be shown by using the fundamental solution of the heat equation after estimating the reaction function f as a non-negative value, compare (Evans 2010). The same applies to a diffusion term ∇·[didiffA(x)∇qi], where A is for all x∈Ω positive definite. Consequently, the non-negativity of the solutions for the cells with only diffusive spreading is assured.

In the case of chemotactic effects, like for the cytotoxic T cells Tc, we show the non-negativity as well. An example of a reaction–diffusion equation for Tc is given byTc,t=FTc(t,x,q1,Tc,q3)+∇·[dTcchemTc∇q3],

analogously for q2. As the reaction function fulfills FTc(t,x,q1,0,q3)≥0, the non-negativity of Tc is proven if the divergence term does not lead to Tc(t,x)<0 for any t>0 and any x∈Ω. The chemotaxis term is the second example in Eq. (18).

Under the assumption of a continuous function, Tc(t,x) is zero before it might become negative. If Tc is negative at a time t and a place x, then the term for the chemotactic effects is negative as well, due to the linear dependency on Tc. So, Tc has non-negative values if the taxis term depends explicitly on Tc, analogously for q2.

The first example in Eq. (18) does not depend on Tc. This term might lead to a negative value of Tc, for example, if Tc(t,x)≡0 for a certain time t and all x∈Ω. Then, depending on the distribution of q3, a negative diffusion Δq3 is possible for a point x. The reaction function can be zero at this point, for example, if q1(t,x)≡0 as well. It is required, that the taxis term depends on the T cellsDTcchem=DTcchem(Tc,∇q3).

Altogether, systems fulfilling the requirements on the reaction functions and on the taxis terms are candidates for well-suited models w.r.t. non-negative solutions. Of course, the requirements are only necessary and not sufficient in all cases.

Boundedness of the solutions

The non-negativity of solutions, which is a lower bound of the solutions, was discussed in Sect. 3.2.1. Now, an upper bound for the solutions is discussed. The effort to show the existence of an upper bound of the solutions depends strongly on the chosen reaction functions. Of course, mechanisms with a negative impact on the variable do not increase an upper bound. If the mechanisms with a positive impact on the variable are already bounded, like some examples in Eqs. (10), (11), and (13), then the solution is bounded. If one mechanism is unbounded, proving an upper bound requires more effort.

From a biological perspective, the existence of upper bounds is a desired property. As the space in the liver is finite, there cannot be infinitely many particles.

We discuss the boundedness of q1, q2, and q3 separately. The existence of upper bounds for Th and Tc follows along the discussion of q2. In some cases, an upper bound of one variable is already required for formulating some other mechanisms.

Again, we start with the boundedness of the virus q1. As the virus needs liver cells for reproduction, the reproduction is limited by the number of liver cells in the domain Ω. Additionally, there can be free virus in the extracellular space, but again, this space is limited as well. Consequently, a growth function of mechanism M1 f1M1(t,x,q1) should fulfill the requirement that35 f1M1(t,x,C1)=0

for a capacity C1>0. Only the first example in Eq. (10) does not fulfill this requirement. This example models an unbounded growth, which is not biologically realistic.

The mechanism M6 depends on the existence of an upper bound for the virus.

Next, we regard the cytokines q3. The reaction function consists of an increase and a natural decay. The growth function f3M3 may depend on the virus, the T (helper) cells, and the cytokines themselves.

We start discussing the cases of f3M3=f3M3(q1). If, as discussed before, the virus is bounded by C1, then the influence of q1 on the increase of q3 is bounded as well. Together with the natural decay, there will be a value q3(t,x)=C3 where the natural decay equalizes the production of cytokines.

For functions f3M3=f3M3(q1,q2,q3), the influence of q3 gives an upper bound through a bounded increase. In the case of f3M3=f3M3(q1,q3), this upper bound is not required due to the interplay with the natural decay.

So far, we did not discuss the boundedness of q2. If q2(t,x)≤C2 for all t>0 and all x∈Ω, then there is an upper bound for q3 as well. This follows by analogous arguments as for the dependency of the increase only on the virus. If q2 is unbounded, the increase in the quantity of cytokines q3 is unbounded as well.

There are two cases of increasing T cells regarding the boundedness, compare Eq. (11). Either, the mechanism for the increase already includes a point-wise upper bound C2 or not. If the mechanism does not include such as boundedness, then the discussion needs more effort, compare (Reisch and Langemann 2022). We refer to this longer discussion of a maximal value and formulate as a requirement that there is a value C~2 such that the time derivative q2,t(t,x) becomes negative for q2(t,x)>C~2.

This requirement is for example fulfilled by the model in Remark 3, see Reisch and Langemann (2022).

Altogether, for gaining a realistic model, all variables should be in L∞(ΩT). If the boundedness of the virus is included in its growth function, the boundedness of the T cells q2 and the cytokines q3 follow in some cases.

For more general models, the regularity and boundedness of solutions for a non-local reaction–diffusion equation is still an active field of research. There is growing knowledge on handling non-local terms in the sense of convolutions, in particular for non-local diffusion, compare (Carrillo et al. 2017). Non-local reaction functions are less frequently studied but biologically relevant, compare (Banerjee et al. 2022).

Longtime behavior

As the solution q=q(t,x) gives the amount of virus, T cells and cytokines at a time t and a space x, the time evolution of the solution can be interpreted as an infection course. The question answered in this section is what types of solutions can be expected and whether these can be predicted a priori. Solutions of reaction–diffusion equations show different behavior like a tendency towards steady states, traveling waves, or blow-ups, compare (Perthame 2015).

In the light of modeling inflammations, requirements for the boundedness of solutions were formulated. These requirements and a bounded domain Ω lead to the boundedness of solutions in any Lp-norm. Consequently, blow-ups are neither desirable for modeling inflammations nor occurring. Solutions of traveling waves are common for modeling atherogenesis, a certain form of inflammation, compare (Ibragimov et al. 2005; Volpert 2014). This type of solution is not suitable for modeling liver infections with a chronic infection course, due to the biological observations of more stationary inflammations, compare (Schiff et al. 2018). The observed spread of T cells can be interpreted as a spatially inhomogeneous stationary solution.

Remark 4

We interpret solutions tending towards a stationary state which is spatially inhomogeneous as chronic liver infections, compare (Kerl et al. 2012). Solutions tending towards zero are interpreted as healing infection courses.

Suitable models for inflammations should show both relevant solution types, depending on the parameters used in the model or on the domain Ω. Analytically, it is desirable to have a priori knowledge about the long-term behavior of the solution. As the tendency towards a spatially homogeneous steady state is a well-analyzed behavior of reaction–diffusion equations, compare for example (Smoller 1994), there is some hope to predict those solutions. The theory of reaction–diffusion equations with solutions tending towards a spatially inhomogeneous steady-state solution is still an open research field.

We discuss two approaches for a priori statements on the occurrence of decaying solutions. The first approach (Smoller 1994, Thm. 14.17), was used in Kerl et al. (2012) in the context of a selected model for inflammations.

The theorem in Smoller (1994) yields for reaction–diffusion equations q,t=F(q)+DΔq. This theorem is not directly applicable to models with chemotactic effects. The statement of the theorem requires some definitions.

Definition 1

(Compare Smoller 1994) Let λ be the first non-zero eigenvalue to the eigenfunctions of the negative Laplacian on Ω with zero flux boundary conditions. Further, let d be the smallest eigenvalue of the diffusion matrix D. Let Σ be the invariant domain of the solutions q, M=maxq∈Σ‖∂F∂q‖Rn×n and σ=λd-M.

Remark 5

(Compare Smoller 1994) The reaction-diffusion system q,t=F(q)+DΔq has an invariant domain of the solutions Σ, F is smooth and D is positive definite. If σ>0 yields, then‖q(t,x)-q¯(t)‖L∞(Ω)≤ce-σt

where q¯(t) is the solution of q¯,t=F(q¯)+g(t) with g decaying exponentially.

Further, if D is diagonal, then q tends in L∞(Ω) towards q¯.

Similar results with slightly different requirements can be found in Perthame (2015).

Next, the requirements of Remark 5 are checked for the models for inflammations. As the structure of the reaction–diffusion equations only considers diffusion, models with chemotactic effects are not covered by the theorem. This is a problem for all models with cytokines. Models with mechanisms like in the first example of Eq. (18) can be written in the form of Remark 5 but those do not fulfill the conditions for non-negative solutions.

Solutions interpreted as chronic infections may occur if σ<0. This is a necessary condition for the occurrence of solutions tending toward a spatially inhomogeneous steady-state distribution. In (Kerl et al. 2012) this is proven for a model similar to the one in Remark 3. In (Reisch and Langemann 2019a) some simplified models with chemotactic effects are analyzed.

For the wide class of possible models in this paper, the theorem of Remark 5 is not well suited and fits only for models without chemotactic effects. Extensions of this theorem or related ones to a wider class of systems are desired.

A different approach for predicting the longtime behavior are entropy methods, compare (Jüngel 2016). As those methods require the definition of an entropy functional depending on the reaction and diffusion terms, a general statement on the use of entropy functionals is not possible. In (Reisch and Langemann 2020) a simplified model for liver inflammations consisting only of one variable was analyzed using entropy methods. The main analytical result for the reduced model is that if the decay terms of the virus overrule the growth of the virus, then only solutions interpreted as healing infection courses occur. This result is not directly transferable to the analysis of a two or more component model because the interplay between the different components dominates the system behavior.

Altogether, only a few conclusions can be drawn from the existing theorems on the longtime behavior of very general reaction–diffusion systems. First, models fulfilling the requirements of Remark 5, that means models without chemotactic effects but with smooth reaction functions with respect to q and upper bounds for all variables, can be analyzed using the theorem from Smoller (1994, Thm. 14.17). The theorem offers two options concerning the two solution types of interest. On the one hand, by regarding as well the ordinary differential equation for the mean q¯, conclusions on the tendency towards a spatial homogeneous solution are possible. On the other hand, like used in Kerl et al. (2012), the contradiction of the theorem can be used, and it can be argued, that in certain cases solutions interpreted as chronic infections are at least possible.

Second, entropy methods turn out to be less effective for general classes of models, especially for models containing explicit space-depending terms, see Reisch and Langemann (2020). If there is an entropy functional for a model, a prediction of a relaxing solution is possible. If the existence of such an entropy functional does not depend on any parameter values, only healing infection courses can be modeled. An open mathematical question therefore is whether entropy-based methods can be used in complex and non-local, space-depending settings.

The next section gives an overview of the requirements of the reaction and taxis functions. Additionally, the requirements are checked for the three presented models.

Summarize of the requirements

The requirements are summarized according to the cell types and the introduced mechanisms. First, the requirements R.1 belonging to the virus are presented.

Remark 6

The dynamics of the virus are described by36 q1,t=f1M1(t,x,q1)+f1M5(t,x,q1,q2)+∇·D1(x,∇q1)

where f1M1, f1M5 and D1 fulfill the requirements (R.1.1) f1,q1M1≥0 for q1>ε>0, compare Eq. (4),

(R.1.2) f1M1(t,x,0)=0, compare Eq. (29),

(R.1.3) f1M1(t,x,C1)=0 for a capacity C1>ε>0, compare Eq. (35)

(R.1.4) f1M5(t,x,0,q2)≥0, compare Eq. (5),

(R.1.5) f1,q2M5≤0, compare Eq. (29),

(R.1.6) if D has the form A∇q1, then A is positive definite.

If the general T cells are divided into T helper cells and cytotoxic T cells, the dependency of f2M5 on q2 is replaced by a dependency on the cytotoxic T cells Tc.

Remark 7

The dynamics of the T cells are described by37 q2,t=f2M2(t,x,q1,q2)+f2M6(t,x,q1,q2)+∇·D2(x,q2,∇q2,∇q3)

which can be adapted for only T helper cells (Th instead of 2 as index) or cytotoxic T cells (Tc instead of 2 as index). In any case, the functions f2M2, f2M6 and D2 fulfill the requirements (R.2.1) f2,q1M2≥0, compare Eq. (6),

(R.2.2) f2M2(t,x,q1,q2)≥0,

(R.2.3) f2,q2M6(t,x,q1,q2)≤0, compare Eq. (6)

(R.2.4) f2M6(t,x,q1,0)≥0, compare Eq. (32),

(R.2.5) there exists a value C2 with q2,t≤0 for all q2>C2,

(R.2.6) if D2 has a function D2ecs(t,x,∇q2) with A∇q2, then A is positive definite,

(R.2.7) if D2 has a function D2chem(q2,∇q3), then the dependency on q2 is obligatory.

As a last quantity, the cytokines are investigated.

Remark 8

The dynamics of the cytokines q3 are described by38 q3,t=f3M3(t,x,q1,q2,q3)+f3nd(t,x,q3)+∇·D3(x,∇q3)

where f3M3, f3nd and D3 fulfill the requirements (R.3.1) f3,q1M3≥0 and f3,ThM3≥0, compare Eq. (7),

(R.3.2) f3M3(t,x,q1,q2,q3)≥0,

(R.3.3) f3nd(t,x,q3)≤0 for all q3≥0,

(R.3.4) f3nd(t,x,0)=0, compare Eq. (34),

(R.3.5) there exists a value C3 with q3,t≤0 for all q3>C3,

(R.3.6) if D3 has the form A∇q3, then A is positive definite.

If the T cells are not split up, the partial derivative in (R.3.1) with respect to Th is replaced by the partial derivative with respect to q2.

The requirements in Remarks 6, 7, 8 are reasonable conditions on the function of a thinkable model for inflammations. Of course, they are neither sufficient nor necessary for formulating a model describing inflammation in a general context.

For the first application, the requirements are checked for the models 1,2 and 3.

Remark 9

Model 1 (Remark 1) fulfills all requirements (R.1), (R.2), (R.3).

Remark 10

Model 2 (Remark 2) fulfills all requirements (R.1), (R.2) and (R.3) under the assumption that q1≤C1. Due to the requirements, the boundedness of the virus is preserved by the dynamics if the initial conditions are locally bounded by C1. The boundedness of the T cells q2 is a consequence of the observations on the two-component model in Reisch and Langemann (2022).

Remark 11

Model 3 in Remark 3 fulfills the requirements (R.1) and (R.2). The boundedness of the T cells q2 was shown in Reisch and Langemann (2022).

In total, all three models fulfill the requirements, even though, in some cases, the proof of the boundedness is more complicated and a priori predictions are not always possible.

Numerical results

Simulations provide a visual understanding of models and their solutions. In this section, simulations for the models 1–3 are presented.

The simulations use a semi-discretization of the space coordinates for Ω=(0,1)×(0,1)⊂R2. Afterward, standard solvers for ordinary differential equations are used. As the motion terms include Laplacians, the resulting ordinary differential equations of the space discretization may become stiff due to the growing quotient of the smallest and the largest eigenvalue of the discretized Laplacian. The number of discretization points is fixed by Δx=0.05 and, in case of long calculation time for solving the ordinary differential equations, a solver for stiff equations is used.

This section focuses on two points. First, the two solution types interpreted as chronic and as healing infection courses are presented. Second, there is a comparison of the three models using two to four different quantities. A leading question is whether even a model covering only two of four cell types can provide any insight.

The initial conditions are chosen as q1(0,x)≡1, q2(0,x)=Tc(0,x)=Th(0,x)≡0 and q3(0,x)≡0.1 in all simulations. The function χΘ(x) is the normalized characteristic function for Θ. Most of the parameters are fixed for all simulations, see Table 1.Table 1 Parameter values for all three models, which are chosen constant if the mechanism is included

a1	C1	ε	κ	d1ctc	a2,h	CTh	a6	dThdiff	CTc	a3	and	d3diff	
1	1	0.05	0.01	0.6	2	8	0.2	0.9	15	0.8	0.6	0.5	

Model 1 shows solutions with a tendency towards zero and solutions with a tendency towards a spatial inhomogeneous steady state, depending on one parameter. In this case, the parameter dTcdiff changes the behavior, see Fig. 2 and Table 2 for the parameter values used. The spread of the virus is modeled by diffusive cell-to-cell transmission. Extracellular movement through a blood flow is not considered, d1ecs=0.Fig. 2 Model 1 shows solutions interpreted as healing infection courses (upper row) and as chronic infections (lower row). Different time steps are in different shadings from t0=0 in light over t1=16 to t2=40 in dark

Table 2 Parameter values for the changed parameters in bold

	Model 1	Model 2	Model 3	
	a5	a2,c	dTcchem	a5	a2	d2chem	a5	a2	
Healing	2	2	1	1	2	1	0.5	2	
Chronic	2	2	8	0.5	2	1	0.5	0.7	

The simulations of model 1 show that there is a remarkable difference in the modeled spreading of T helper cells Th and cytotoxic T cells Tc. While the spread of T helper cells is modeled by diffusion in model 1, compare Remark 1, the motion of cytotoxic T cells is modeled by chemotaxis. Consequently, the cytotoxic T cells Tc in Fig. 2g spread out less and the immune reaction is less effective in reducing the viral load. Figure 2c, g show the inflow region Θ of the immune response, which is located in the lower-right corner of the spatial domain. Due to the chemotaxis, the simulation of the cytotoxic T cells Tc show spiky behavior in the region Θ, where the increase of Tc depending on the virus u happens. In the chronic case, the number of T helper cells Th reduces slightly after an active phase but remains high for all time. In contrast, the number of T helper cells Th decays to zero in the healing infection course once the virus is reduced.

In model 2, there is only one type of T cells spreading out by diffusion and following the gradient of the cytokines. The parameter under variation in this case is a5, see Table 2, which describes the effectiveness of the T cells. Figure 3 shows a solution interpreted as chronic infection course. The virus remains in the whole liver with a spatial inhomogeneous distribution. The cytokines remain at a high level everywhere.Fig. 3 Model 2 shows solutions interpreted as chronic infection courses. Different time steps are in different shadings from t0=0 in light over t1=32 to t2=80 in dark. The portal field Θ is in the right corner in (b)

The changed parameter for model 3 is a2, see Table 2, which regulates the inflow of T cells through the portal field in dependency on the total virus at a time t. The results of the simulation of model 3 are comparable to the simulation results of models 1 and 2. A comparison of all three models for a chronic course is in Fig. 4. As the T cells q2 are a sum of the T helper cells Th and the cytotoxic T cells Tc, Fig. 4 compares the L1 norms of q1 in (b) of q2 and Th+Tc.Fig. 4 L1 norms of the solutions of model 1, 2 and 3 for the parameter sets in Table 2 interpreted as chronic infection course. a Total virus, b total T cells

The amount of virus tends in all three models to comparable values due to the small changes of parameters. The total amount of T cells is comparable for models 2 and 3 but different for model 1 due to the differentiation between T helper cells and cytotoxic T cells with different maximal values. The curves of models 2 and 3 show small oscillations in the first time but a damping towards a constant value later on.

Clinical impact

The mathematical model for liver inflammation developed and described here may have significant clinical implications, as it provides valuable insights into the underlying dynamics and mechanisms of liver inflammation caused by viruses. Some potential clinical improvements include a more profound understanding of the disease and its progression, supporting clinical decisions through early detection and diagnosis, and optimizing treatment. All of this could contribute to personalized diagnosis and treatment of viral hepatitis in the future. More specifically, this mathematical model can help to detect subtle changes in influencing variables or physiological parameters associated with liver inflammation at an early stage. This early detection allows for rapid intervention and treatment that can potentially prevent the progression of liver inflammation to cirrhosis. Such predictive models can help physicians use their resources more efficiently by identifying patients who are likely to require more intensive monitoring or intervention. This is of paramount importance in healthcare systems with limited resources. Mathematical models provide a platform for simulating the dynamic processes involved in liver inflammation. Understanding disease progression at a mechanistic level can contribute to the development of new therapeutic targets and strategies, or prevent inflammation from becoming chronic. More and more patients are demanding shared decision-making, where the visualization and explanation derived from the mathematical model can be used to educate patients about their disease. Better patient understanding can improve treatment adherence and lifestyle change, which has a positive impact on long-term outcomes. In summary, a mathematical model for liver inflammation can provide clinicians with valuable tools to improve diagnosis, treatment planning, and patient outcomes. Its integration into clinical practice has the potential to usher in a new era of precision medicine for liver disease.

Conclusions

A deductive modeling approach was presented for the application of liver infections leading to inflammations. This life science application has unknown mechanisms leading to chronic infection courses, for which only qualitative data is available. Building up a model family allows gaining a more profound understanding of the involved mechanisms. The model family consists of reaction–diffusion equations, and the reaction or taxis functions can be chosen from a class of feasible functions. The feasible classes fulfill different biologically motivated properties. Analytical investigations restricted the function classes further.

In the case of inflammation modeling, three models were chosen out of the model family and all three models reproduce both different infection courses depending on some parameter values. Consequently, already the smallest of the three models can cover basic observations even if less information is used than in the more complex models of the model family. Depending on the modeling purpose, a suitable model from the model family can be chosen.

Concerning the analysis of reaction–diffusion equations with space-dependent non-local terms, there are some open questions: The parameter value of the space-dependent non-local term, here named a2, is a bifurcation parameter for the steady state, compare Table 2. An open question is the prediction of the change of the solution behavior, depending on the parameter and the exact non-local function. Analytical approaches in the literature (Van Gorder 2021) are based on the linearization of the equations, Fourier-approaches for the solutions and truncations of the resulting infinite-dimensional linear system for the Fourier coefficients. One challenge is the treatment of the coupling of the coefficients due to the space-dependent and non-local term. This strong coupling complicates any truncation of the infinite dimensional system and the study of the stability of single Fourier modes.

The stability of stationary states, and in particular of spatially heterogeneous steady states, is an open question that might be answered by using computer-assisted proofs for studying the convergence of the dynamical system towards those states. The heterogeneous steady states are interpreted as chronic inflammations. The development of those permanently inflammatory states is medically relevant, therefore the analytical question of convergence towards such states is not only challenging but interesting in application.

A more profound understanding of the longtime behavior of the solutions allows predicting a priori the influence of parameter changes. In the clinical context, medication leads to changes of the parameters. The combination of mathematical models and analytical predictions may help to personalize the medical treatment.

Funding

Open Access funding enabled and organized by Projekt DEAL. HM Tautenhahn disclosed receipt of the following financial support for the authorship of this article: This work was supported by: Federal Ministry of Education and Research; ATLAS (031L0304C) HM Tautenhahn German Research Foundation SIMLIVA (465194077) HM Tautenhahn German Research Foundation (428832822) HM Tautenhahn This research did not contain any studies involving animal or human participants, nor did it take place in any private or protected areas. No specific permissions were required for corresponding locations.

Data availibility statement

The paper is not using experimental data sets.

Declarations

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

Allee WC Principles of animal ecology 1949 Philadelphia Saunders CO
Allee WC (1949) Principles of animal ecology. Saunders CO, Philadelphia
Aston P A new model for the dynamics of hepatitis C infection: derivation, analysis and implications Viruses 2018 10 4 195 10.3390/v10040195 29652855
Aston P (2018) A new model for the dynamics of hepatitis C infection: derivation, analysis and implications. Viruses 10(4):195. 10.3390/v1004019529652855 10.3390/v10040195
Banerjee M Kuznetsov M Udovenko O Nonlocal reaction–diffusion equations in biomedical applications Acta Biotheoretica 2022 70 2 12 10.1007/s10441-022-09436-4 35298702
Banerjee M, Kuznetsov M, Udovenko O et al (2022) Nonlocal reaction–diffusion equations in biomedical applications. Acta Biotheoretica 70(2):12. 10.1007/s10441-022-09436-435298702 10.1007/s10441-022-09436-4
Bowen DG Walker CM Adaptive immune responses in acute and chronic hepatitis C virus infection Nature 2005 436 7053 946 952 10.1038/nature04079 16107834
Bowen DG, Walker CM (2005) Adaptive immune responses in acute and chronic hepatitis C virus infection. Nature 436(7053):946–952. 10.1038/nature0407916107834 10.1038/nature04079
Carrillo JA, Del Pino M, Figalli A et al (2017) Nonlocal and Nonlinear Diffusions and Interactions: New Methods and Directions, Lecture Notes in Mathematics. Springer, Cham. 10.1007/978-3-319-61494-6
Evans LC Partial Differential Equations, Graduate studies in mathematics 2010 2 Providence, RI American Mathematical Society
Evans LC (2010) Partial Differential Equations, Graduate studies in mathematics, vol 19, 2nd edn. American Mathematical Society, Providence, RI
Goyal A Murray JM Modelling the impact of cell-to-cell transmission in hepatitis B virus PLoS ONE 2016 11 8 e0161978 10.1371/journal.pone.0161978 27560827
Goyal A, Murray JM (2016) Modelling the impact of cell-to-cell transmission in hepatitis B virus. PLoS ONE 11(8):e0161978. 10.1371/journal.pone.016197827560827 10.1371/journal.pone.0161978
Ibragimov AI McNeal CJ Ritter LR A mathematical model of atherogenesis as an inflammatory response Math Med Biol 2005 22 4 305 333 10.1093/imammb/dqi011 16162594
Ibragimov AI, McNeal CJ, Ritter LR et al (2005) A mathematical model of atherogenesis as an inflammatory response. Math Med Biol 22(4):305–333. 10.1093/imammb/dqi01116162594 10.1093/imammb/dqi011
Jüngel A Entropy methods for diffusive partial differential equations 2016 Cham Springer
Jüngel A (2016) Entropy methods for diffusive partial differential equations. Springer, Cham
Kanel GC Pathology of liver diseases 2017 Hoboken, NJ Wiley
Kanel GC (2017) Pathology of liver diseases. Wiley, Hoboken, NJ
Kazankov K Jørgensen SMD Thomsen KL The role of macrophages in nonalcoholic fatty liver disease and nonalcoholic steatohepatitis Nat Rev Gastroenterol Hepatol 2019 16 3 145 159 10.1038/s41575-018-0082-x 30482910
Kazankov K, Jørgensen SMD, Thomsen KL et al (2019) The role of macrophages in nonalcoholic fatty liver disease and nonalcoholic steatohepatitis. Nat Rev Gastroenterol Hepatol 16(3):145–159. 10.1038/s41575-018-0082-x30482910 10.1038/s41575-018-0082-x
Kerl HJ Langemann D Vollrath A Reaction diffusion equations and the chronification of liver infections Math Comput Simul 2012 82 11 2145 2156 10.1016/j.matcom.2012.04.011
Kerl HJ, Langemann D, Vollrath A (2012) Reaction diffusion equations and the chronification of liver infections. Math Comput Simul 82(11):2145–2156. 10.1016/j.matcom.2012.04.01110.1016/j.matcom.2012.04.011
Murray JD Mathematical Biology, Interdisciplinary applied mathematics. Spatial Models and Biomedical Applications 2002 3 New York Springer
Murray JD (2002) Mathematical Biology, Interdisciplinary applied mathematics. Spatial Models and Biomedical Applications, vol II, 3rd edn. Springer, New York
Nangue A Rendall AD Tcheugam BK Analysis of an initial value problem for an extracellular and intracellular model of hepatitis C virus infection Int J Biomath 2022 10.1142/S1793524522500413
Nangue A, Rendall AD, Tcheugam BK et al (2022) Analysis of an initial value problem for an extracellular and intracellular model of hepatitis C virus infection. Int J Biomath. 10.1142/S179352452250041310.1142/S1793524522500413
Perthame B Parabolic equations in biology: growth, reaction, movement and diffusion. Lecture notes on mathematical modelling in the life sciences 2015 Cham Springer
Perthame B (2015) Parabolic equations in biology: growth, reaction, movement and diffusion. Lecture notes on mathematical modelling in the life sciences. Springer, Cham
Reisch C Langemann D Chemotactic effects in reaction–diffusion equations for inflammation J Biol Phys 2019 45 3 253 273 10.1007/s10867-019-09527-3 31309352
Reisch C, Langemann D (2019) Chemotactic effects in reaction–diffusion equations for inflammation. J Biol Phys 45(3):253–273. 10.1007/s10867-019-09527-331309352 10.1007/s10867-019-09527-3
Reisch C Langemann D Modeling the chronification tendency of liver infections as evolutionary advantage Bull Math Biol 2019 81 11 4743 4760 10.1007/s11538-019-00596-y 30891676
Reisch C, Langemann D (2019) Modeling the chronification tendency of liver infections as evolutionary advantage. Bull Math Biol 81(11):4743–4760. 10.1007/s11538-019-00596-y30891676 10.1007/s11538-019-00596-y
Reisch C Langemann D Entropy functionals for finding requirements in hierarchical reaction–diffusion models for inflammations Math Meth Appl Sci 2020 43 17 10098 10114 10.1002/mma.6682
Reisch C, Langemann D (2020) Entropy functionals for finding requirements in hierarchical reaction–diffusion models for inflammations. Math Meth Appl Sci 43(17):10098–10114. 10.1002/mma.668210.1002/mma.6682
Reisch C Langemann D Longterm existence of solutions of a reaction diffusion system with non-local terms modeling an immune response—an interpretation-orientated proof Partial Differ Equ Appl Math 2022 6 100446 10.1016/j.padiff.2022.100446
Reisch C, Langemann D (2022) Longterm existence of solutions of a reaction diffusion system with non-local terms modeling an immune response—an interpretation-orientated proof. Partial Differ Equ Appl Math 6:100446. 10.1016/j.padiff.2022.10044610.1016/j.padiff.2022.100446
Rezounenko A Viral infection model with diffusion and state-dependent delay: stability of classical solutions Discrete Contin Dyn Syst B 2018 23 3 1091 1105 10.3934/dcdsb.2018143
Rezounenko A (2018) Viral infection model with diffusion and state-dependent delay: stability of classical solutions. Discrete Contin Dyn Syst B 23(3):1091–1105. 10.3934/dcdsb.201814310.3934/dcdsb.2018143
Schiff ER, Maddrey WC, Reddy KR (eds) (2018) Schiff’s diseases of the liver, twelfth edn. Wiley, Chichester
Schnakenberg J Simple chemical reaction systems with limit cycle behaviour J Theor Biol 1979 81 3 389 400 10.1016/0022-5193(79)90042-0 537379
Schnakenberg J (1979) Simple chemical reaction systems with limit cycle behaviour. J Theor Biol 81(3):389–400. 10.1016/0022-5193(79)90042-0537379 10.1016/0022-5193(79)90042-0
Smoller J Shock waves and reaction–diffusion equations 1994 2 New York Springer
Smoller J (1994) Shock waves and reaction–diffusion equations, 2nd edn. Springer, New York
Tadmon C Foko S Rendall AD Global stability analysis of a delay cell-population model of hepatitis B infection with humoral immune response Dyn Syst 2021 36 4 537 559 10.1080/14689367.2021.1940868
Tadmon C, Foko S, Rendall AD (2021) Global stability analysis of a delay cell-population model of hepatitis B infection with humoral immune response. Dyn Syst 36(4):537–559. 10.1080/14689367.2021.194086810.1080/14689367.2021.1940868
Thomas E Liang TJ Experimental models of hepatitis B and C—new insights and progress Nat Rev Gastroenterol Hepatol 2016 13 6 362 374 10.1038/nrgastro.2016.37 27075261
Thomas E, Liang TJ (2016) Experimental models of hepatitis B and C—new insights and progress. Nat Rev Gastroenterol Hepatol 13(6):362–374. 10.1038/nrgastro.2016.3727075261 10.1038/nrgastro.2016.37
Turing A The chemical basis of morphogenesis Philos Trans R Soc Lond B Biol Sci 1952 237 641 37 72 10.1098/rstb.1952.0012
Turing A (1952) The chemical basis of morphogenesis. Philos Trans R Soc Lond B Biol Sci 237(641):37–72. 10.1098/rstb.1952.001210.1098/rstb.1952.0012
Van Gorder RA Pattern formation from spatially heterogeneous reaction–diffusion systems Philos Trans R Soc A: Math Phys Eng Sci 2021 379 2213 20210001 10.1098/rsta.2021.0001
Van Gorder RA (2021) Pattern formation from spatially heterogeneous reaction–diffusion systems. Philos Trans R Soc A: Math Phys Eng Sci 379(2213):20210001. 10.1098/rsta.2021.000110.1098/rsta.2021.0001
Volpert V Elliptic partial differential equations, vol 2: reaction–diffusion-equations, monographs in mathematics 2014 Basel Springer
Volpert V (2014) Elliptic partial differential equations, vol 2: reaction–diffusion-equations, monographs in mathematics, vol 104. Springer, Basel
