A Mathematical Model for Transport in Poroelastic Materials with Variable Volume: Derivation, Lie Symmetry Analysis and Examples—Part 2

Abstract

The model for perfused tissue undergoing deformation taking into account the local exchange between tissue and blood and lymphatic systems is presented. The Lie symmetry analysis in order to identify its symmetry properties is applied. Several families of steady-state solutions in closed formulae are derived. An analysis of the impact of the parameter values and boundary conditions on the distribution of hydrostatic pressure, osmotic agent concentration and deformation of perfused tissue is provided applying the solutions obtained in examples describing real-world processes.

1. Introduction

For any living systems, fluid and solute transport are crucial for their functioning. Despite complexity of the regulatory processes, the natural tendency of living systems to maintain optimal condition of functioning brings them to the homeostatic state (local steady-state). Any perturbation of this state, by external or internal forces, triggers regulatory processes to get a new steady-state. In this study, we focus at the transport processes occurring in the biological tissue, such as the muscle tissue or tumour tissue [1,2,3]. Due to the complexity of the processes and nonlinear interactions between system components, the mathematical description of the transport processes has to be used for better understanding of the underlying physiology and pathology of the observed phenomenae.

The local tissue properties are not constant but they are changing dynamically altering fluid and solute transport through the tissue [4,5,6,7,8]. Changes in the tissue hydration are closely link to the changes of the interstitial pressure by modifying some tissue transport parameters, such as tissue hydraulic conductance, solute diffusivity in the tissue or lymphatic absorption, as shown on the ground of local tissue physiology in various experiments [9,10]. Moreover, in the case of insufficient regulatory mechanisms, such changes may also lead to the volume and shape modifications, that are rarely taken into account in the mathematical models [4,6,7]. Such simplified approach can be used in the case of small perturbation but is not valid in the case of larger perturbation caused, for example, by the substantial change of the external forces.

The modeling of combined transport and elastic deformation phenomenae in physiology are rather rare in general and often focus at the tissue that is not perfused by blood [11,12,13,14], except for a series of papers on solid tumors [4,15,16]. Even then the deformation is often assumed small and neglected [4,6,16]. All these problems get more sophisticated if the soft tissue perfused by blood and lymphatic vessels is considered because of the additional exchange of fluid and solute between the tissue and blood and lymph [3,4,15,16,17]; in particular, the inflow of water to the tissue from the peritoneal cavity due to high hydrostatic pressure of dialysis fluid and from the capillaries due to the increased osmotic pressure of interstitial fluid (tissue fluid phase) [7,18,19,20]. The overhydration of tissue may result in deformation of the tissue layer and its increased width. In general, the tissue may also shrink under, for example, decreased external hydrostatic pressure due to the removal of fluid or external osmotic pressure without the increase in hydrostatic pressure.

Our previous analysis considered nonperfused poroelastic material (PEM). The linear theory of poroelasticity was applied in order to model its volume/shape deformation due to external hydrostatic and/or osmotic pressure [21]. In this study, we extend previously developed model focusing on the biological tissue perfused by blood and drained by lymphatic system [10,17,22]. An important feature of our approach is the role of osmotic pressure that is not considered in nonbiological applications of elastic theory, as in geology for the description of soil or rock perfused by water [14,23,24]. Osmotic pressure is crucial for the distribution and shifts of water in the body [1], and is the main (thermodynamic) force for the removal of excess of water during peritoneal dialysis of renal failure patients [7,18]. The inclusion of osmotic pressure into the theory leads to increase complexity of the model by adding nonlinear terms into transport equations. In the study we applied the concept of tissue deformation – one-dimensional poroelastic theory with the extended Terzaghi stress tensor that takes into account osmotic effect [23]. Moreover, the generalized nonlinear description of stress tensor with a quadratic term is considered.

In Section 2, the model for perfused tissue undergoing deformation taking into account the local exchange between tissue and blood and lymphatic systems is presented. The model consists of five nonlinear PDEs, which should be supplied by the relevant boundary conditions. Since the model is very complex, we used the classical Lie method [25,26,27], which is one of the most powerful tools for constructing exact solutions of nonlinear PDEs, in order to identify its Lie symmetry properties. Although the Lie symmetries of the model are rather trivial in the general case, some special cases were identified when additional nontrivial symmetries appear.

In Section 3, a nonlinear system of ordinary differential equations (ODEs) corresponding to the steady-state of the model is under study. Since the ODE system is still nonintegrable we demonstrate that its nontrivial exact solutions can be constructed under correctly-specified restrictions on its parameters. In particular, exact solutions involving several arbitrary constants are derived.

In Section 4, the formulae derived in previous sections have been applied to study deformation of a perfused tissue. The new equilibrium of hydrostatic pressure profile, solute concentration and tissue deformation caused by external force, which leads to the defined stretching, has been simulated numerically. We show how the elastic modulus and nonlinearity of the stress tensor influence deformation of the porous tissue. Moreover, impact of external pressure on the pressure profiles and deformation profiles has been provided.

2. Model for Perfused Tissue Undergoing Deformation and Its Lie Symmetry

Let us consider the steady-states equations for fluid and solute transport though the tissue. The fluid flux is composed of the volumetric flux relative to the matrix, driven by the extended Darcy’s law, and fluid flow between tissue and circulatory system, $q_V$. Let us assume that local exchange occurs only with blood and lymphatic systems, so that fluid inflow into the tissue is

$$q_V=q_C-q_L,$$

where $q_L$ denotes density of volumetric fluid absorption and fluid flows across blood capillary wall (formally density of the volumetric fluid flow) is driven according to the Starling’s law as follows

$$q_C=a{l}_P(P_B-p-{\sigma }_{C}RT(C_B-c)),$$

where $p(t,x)$ and $c(t,x)$ are tissue (interstitial) hydrostatic pressure and solute concentration, respectively, being functions of time t and position x, and $a{l}_P,{\sigma }_C,P_B,RT,C_B$ are fixed constants, see

Table 1. Description of the symbols.

Symbol	Description
c	solute concentration in the tissue
p	interstitial hydrostatic pressure in the tissue
u	deformation vector
$\rho$	mass density
${\rho }_F$	mass density of fluid phase F
${\rho }_M$	mass density of matrix phase M
${\theta }_F$	fractional volume of fluid in fluid phase F
${\theta }_M$	fractional volume of fluid in solid phase M
$q_L$	density of volumetric fluid absorption by lymphatics
$q_0$	rate of lymphatic absorption for p = 0
$q_1$	sensitivity of lymphatic absorption to increase in p
$q_C$	density of volumetric fluid flow across blood capillary wall
$q_V$	net fluid flow into the tissue from internal sources/sinks
$q_S$	net solute flow into the tissue from internal sources/sinks
$C_B$	solute concentration in blood
$C_M$	mean solute concentration across blood capillary wall
$P_B$	blood hydrostatic pressure
$F_1,F_2$	weighting factors for concentration across blood capillary wall
$a{l}_P$	hydraulic conductance for blood capillary wall
$a{p}_D$	solute diffusive permeability of blood capillary wall
${\sigma }_C$	solute reflection coefficient for blood capillary wall
$\sigma$	solute reflection coefficient for tissue
$S=1-\sigma$	sieving coefficient of solute in tissue
$S_C=1-{\sigma }_C$	solute sieving coefficient for blood capillary wall
$RT$	gas constant times temperature
$\tau$	Terzaghi effective stress tensor
$\lambda \ast =\lambda +2\mu$	elastic modulus with Lame constants $\lambda$ and $\mu$
$\kappa$	parameter for nonliear component of stress tensor
k	hydraulic conductivity of tissue
D	solute diffusivity in the tissue

. Moreover, we assume that lymphatic absorption does not have to be constant but may change proportionally to the changes of interstitial pressure [10,22]:

$$q_L=q_0+q_{1}p,$$

where $q_0$ and $q_1$ are fixed constants and nonnegative in biological applications.

The solute flux through the porous medium consists of the solute flux relative to the medium, that can be both diffusive or convective, and the solute exchange between tissue and circulatory systems, $q_S$. We assume that the rate of net solute inflow into the tissue is due to local exchange across blood capillary wall, with solute transport that can be both diffusive and convective, decreased by local solute absorption by lymphatics with fluid bulk flow, namely

$$q_S=a{p}_D(C_B-c)+S_C{C}_{M}a{l}_P\left(P_B-p-{\sigma }_{C}RT(C_B-c)\right)-q_{L}c,$$

where $a{p}_D$ and $S_C$ are constant related to the transport properties of the blood capillary, whereas $C_M$ denotes mean concentration across the blood capillary wall and typically linearly depends on c.

The deformation of the matrix, $u(t,x)$, is directly related to the stress tensor, that in general can be of nonlinear form [23]. In this study a quadratic form of Terzaghi tensor, being the simplest approximation of this nonlinear dependence, has been considered as follows:

$$\tau =-(p-\sigma RTc)+{\lambda }^{\ast }\frac{\partial u}{\partial x}+\kappa {\left(\frac{\partial u}{\partial x}\right)}^2,$$

where $\sigma$ denotes sieving coefficient for tissue, ${\lambda }^{\ast }=\lambda +2\mu$ is the Lame coefficient [6,21] and the parameter $\kappa$ describes the contribution of squared dilatation to the stress tensor. Obviously, setting $\kappa =0$ one obtains the standard, linear expression for $\tau$.

Having the formulae for fluxes $q_V,q_C$, and $q_S$ presented above, we have constructed the model consisting of five nonlinear PDEs in a quite similar way as it was done in [21] (we remind that internal sources/sinks were not taken into account in that paper). The basic equations of the model read as

$$2u_{tx}=k(p_{xx}-\sigma RT{c}_{xx})+q_V,$$

(1)

$$\rho u_{tt}+{\rho }_t{u}_t+\rho u_t{u}_{tx}={\lambda }^{\ast }{u}_{xx}+2\kappa u_x{u}_{xx}-(p_x-\sigma RT{c}_x),$$

(2)

$${\rho }_t+{\rho }_x{u}_t=k({\rho }_F^0-\rho )(p_{xx}-\sigma RT{c}_{xx})+({\rho }_F^0-\rho )q_V,$$

(3)

$${{\theta }_F}_t+{{\theta }_F}_x{u}_t=k\left(1-{\theta }_F\right)(p_{xx}-\sigma RT{c}_{xx})+(1-{\theta }_F)q_V,$$

(4)

$${\left(c{\theta }_F\right)}_t+{\left(c{\theta }_F\right)}_x{u}_t+2c{\theta }_F{u}_{tx}=D{c}_{xx}+kS{\left(c(p_x-\sigma RT{c}_x)\right)}_x+q_S,$$

(5)

where the lower subscripts t and x denote differentiation with respect to these variables. Although $q_V$ and $q_L$ in system (1)–(5) can be of more complicated forms, however we assume in what follows that

$$q_V=a{l}_P(P_B-p-{\sigma }_{C}RT(C_B-c))-q_L,$$

(6)

$$q_S=a{p}_D(C_B-c)+S_C{C}_M(q_V+q_L)-q_{L}c,$$

(7)

$$q_L=q_0+q_{1}p,C_M=F_1{C}_B+F_{2}c(t,x),F_1+F_2=1,$$

(8)

where $q_0$, $q_1$, $F_1$ and $F_2$ are nonnegative constants.

The above equations were derived under the natural restrictions ${\theta }_F+{\theta }_M=1$ and $\rho ={\rho }_F^0{\theta }_F+{\rho }_M{\theta }_M$, which hold for two-phase poroelastic material (tissue). We also assume that the fluid is incompressible therefore ${\rho }_F^0$ is a positive constant. All the notations arising in the above formulae and in those below are explained in Table 1.

Obviously, the nonlinear PDE system (1)–(5) is an essential generalization of that derived in paper [21] (see Equations (21)–(25) therein). Now we present a theorem presenting Lie’s symmetry properties of (1)–(5). One contains Theorem 1 from [21] as a particular result (see item (v) below).

Theorem 1.

The principal algebra of invariance of system (1)–(5) is a three-dimensional Lie algebra generated by the Lie symmetry operators

$$X_1={\partial }_t,X_2={\partial }_x,X_3={\partial }_u.$$

We remind the reader that the notion of the principal algebra of invariance (see, e.g., [27]) means that a given PDE (system of PDEs) admits this algebra for arbitrary set of parameters arising in the PDE (system of PDEs) in question. Moreover, this algebra of invariance cannot be extended by any other Lie symmetry without a restriction on parameters of the given PDE (system of PDEs).

Theorem 2.

The Lie algebra of invariance <$X_1,X_2,X_3$> of system (1)–(5) depending on the parameter κ and the functions $q_V,q_L$ and $q_S$ from (6)–(8) admits the following extensions:

(i)

if $\kappa =0$ then the additional operator is $X_4=x{\partial }_u$;

(ii)

if $\kappa \ne 0$, $q_L=-q_V=q_0$ and $q_S+q_{0}c$ is nonzero then the additional operator is $X_5=g\left(t\right){\partial }_p$;

(iii)

if $\kappa =0$, $q_L=-q_V=q_0$ and $q_S+q_{0}c$ is nonzero then the additional operators are $X_4=x{\partial }_u$ and $X_5=g\left(t\right){\partial }_p$;

(iv)

if $\kappa \ne 0$, $q_L=-q_V=q_0$ and $q_S=-q_{0}c$ then the additional operators are $X_5=g\left(t\right){\partial }_p$ and $X_6=c({\partial }_c+\sigma RT{\partial }_p)$;

(v)

if $\kappa =0$, $q_L=-q_V=q_0$ and $q_S=-q_{0}c$ then the additional operators are $X_4=x{\partial }_u$, $X_5=g\left(t\right){\partial }_p$ and $X_6=c({\partial }_c+\sigma RT{\partial }_p).$

Here $g\left(t\right)$ is an arbitrary smooth function and the notations ${\partial }_z\equiv \frac{\partial }{\partial z}(z=t,x,u,c,p)$ are used.

Sketch of the proof of Theorems 1 and 2. 

It is based on the infinitesimal criteria of invariance, which was formulated by Sophus Lie in his pioneering works published in the end of the 19th century. In the case of a system of PDEs of arbitrary order, this criteria can be found, e.g., in [25] (see Section 1.2.5 therein). Nowadays the criteria is included in several computer algebra packages such as Maple and Mathematica and the result can be derived without long calculations by hands (provided the problem of Lie symmetry classification does not arise). Since the functions $q_L,q_V$ and $q_S$ are explicitly specified in (6)–(8) the problem of Lie symmetry classification reduces to examination of several cases when the parameters arising in (6)–(8) vanish. We used Maple 18 for these purposes and the result is presented in the above theorems. □

It should be noted that the Lie symmetries $X_1,X_2$ and $X_3$ are rather obvious and can be noted without application of the Lie algorithm to system (1)–(5). However, the Lie symmetries $X_4,X_5$ and $X_6$ arising in special cases are rather hidden hence Theorem 2 presents nontrivial results.

Having the Lie symmetry operators defined explicitly, one has several possibilities to reduce system (1)–(5) to systems of ODEs in order to find exact solutions with different structures. In particular, the principal algebra of invariance allows us to construct the linear combination

$$Y={\alpha }_0{\partial }_t+\alpha {\partial }_x+{\alpha }_1{\partial }_u,$$

where ${\alpha }_0,\alpha$ and ${\alpha }_1$ are arbitrary parameters. There are only two essentially different cases: ${\alpha }_0\ne 0$ and ${\alpha }_0=0$. In the first case, we may set ${\alpha }_0=1$, so that the ansatz

$$\begin{array}{c}u=f_1\left(\omega \right)+{\alpha }_{1}t,p=f_2\left(\omega \right),\rho =f_3\left(\omega \right),\hfill \\ {\theta }_F=f_4\left(\omega \right),c=f_5\left(\omega \right),\hfill \end{array}$$

(9)

is obtained. Here $\omega =x-\alpha t$ and $f_i(i=1,\cdots ,6)$ are to-be-determined functions.

Substituting ansatz (9) into system (1)–(5), the system of ODEs

$$\begin{array}{c}k(f_2^{″}-\sigma RT{f}_5^{″})+2\alpha f_1^{″}+q_V=0,\hfill \\ {\lambda }^{\ast }{f}_1^{″}+2\kappa f_1^{\prime }{f}_1^{″}-(f_2^{\prime }-\sigma RT{f}_5^{\prime })+\alpha f_3{f}_1^{″}\left({\alpha }_1-\alpha -\alpha f_1^{\prime }\right)+\alpha ({\alpha }_1-\alpha f_1^{\prime })f_3^{\prime }=0,\hfill \\ k({\rho }_F^0-f_3)(f_2^{″}-\sigma RT{f}_5^{″})-f_3^{\prime }\left({\alpha }_1-\alpha -\alpha f_1^{\prime }\right)+({\rho }_F^0-f_3)q_V=0,\hfill \\ k\left(1-f_4\right)(f_2^{″}-\sigma RT{f}_5^{″})-f_4^{\prime }\left({\alpha }_1-\alpha -\alpha f_1^{\prime }\right)+(1-f_4)q_V=0,\hfill \\ D{f}_5^{″}+kS{\left(f_5(f_2^{\prime }-\sigma RT{f}_5^{\prime })\right)}^{\prime }+2\alpha f_4{f}_5{f}_1^{″}-{\left(f_4{f}_5\right)}^{\prime }\left({\alpha }_1-\alpha -\alpha f_1^{\prime }\right)+q_S=0,\hfill \end{array}$$

(10)

to find the functions $f_i$ is obtained. In system (10):

$$\begin{array}{c}{q}_V=a{l}_P\left(P_B-f_2-{\sigma }_{C}RT(C_B-f_5)\right)-q_L,q_L=q_0+q_1{f}_2,\hfill \\ q_S=a{p}_D(C_B-f_5)+a{l}_P{S}_C(F_1{C}_B+F_2{f}_5)(P_B-f_2-{\sigma }_{C}RT(C_B-f_5))-q_L{f}_5.\hfill \end{array}$$

Setting $\alpha ={\alpha }_1=0$, an ODE system for finding steady-state solutions is obtained, which is studied in Section 3. Another particular case $\alpha \ne 0,{\alpha }_1=0$ leads to the ansatz for search plane wave solutions, which are also useful in applications. Assuming ${\alpha }_0=0$, we may set $\alpha =1$ without losing a generality, so that the ansatz

$$\begin{array}{c}u=f_1\left(t\right)+{\alpha }_{1}x,p=f_2\left(t\right),\rho =f_3\left(t\right),\hfill \\ {\theta }_F=f_4\left(t\right),c=f_5\left(t\right),\hfill \end{array}$$

(11)

is derived. However, exact solutions of the form (11) are not useful from the applicability point of view.

The case (i) arising in Theorem 2 corresponds to the model with the linear Terzaghi tensor, which is commonly used in applications. The relevant Lie algebra is four-dimensional. In contrast to the principal algebra, the algebra <$X_1,X_2,X_3,X_4$> is non-Abelian and leads to a larger number of inequivalent reductions of system (1)–(5) to ODEs. All inequivalent reduction can be identified either via straightforward technique (see, e.g., Section 4.2 [27]), or using an optimal system of inequivalent (nonconjugated) one-dimensional subalgebras. The optimal system was constructed in the seminal work [28] (see the fourth case in Table II) and one may identify that only the linear combination

$$Z={\alpha }_0{\partial }_t+\alpha {\partial }_x+{\alpha }_{1}x{\partial }_u,{\alpha }_1\ne 0$$

leads to new ansatzes comparing with those presented above. So, two essentially different cases again occur: ${\alpha }_0\ne 0$ and ${\alpha }_0=0$.

In the case ${\alpha }_0\ne 0\Rightarrow {\alpha }_0=1$, we derive the ansatz

$$\begin{array}{c}u=f_1\left(\omega \right)+\frac{{\alpha }_1}{2\alpha }{x}^2,p=f_2\left(\omega \right),\rho =f_3\left(\omega \right),\hfill \\ {\theta }_F=f_4\left(\omega \right),c=f_5\left(\omega \right),\omega =x-\alpha t,\hfill \end{array}$$

if $\alpha \ne 0$ and

$$\begin{array}{c}u=f_1\left(x\right)+{\alpha }_{1}tx,p=f_2\left(x\right),\rho =f_3\left(x\right),\hfill \\ {\theta }_F=f_4\left(x\right),c=f_5\left(x\right),\hfill \end{array}$$

(12)

if $\alpha =0$.

The corresponding systems of the reduction equations are

$$\begin{array}{c}k(f_2^{″}-\sigma RT{f}_5^{″})+2\alpha f_1^{″}+q_V=0,\hfill \\ {\lambda }^{\ast }\left(f_1^{″}+\frac{{\alpha }_1}{\alpha }\right)-(f_2^{\prime }-\sigma RT{f}_5^{\prime })-{\alpha }^2\left(f_3{f}_1^{″}+f_1^{\prime }{f}_3^{\prime }+f_3{f}_1^{\prime }{f}_1^{″}\right)=0,\hfill \\ k({\rho }_F^0-f_3)(f_2^{″}-\sigma RT{f}_5^{″})+\alpha f_3^{\prime }\left(1+f_1^{\prime }\right)+({\rho }_F^0-f_3)q_V=0,\hfill \\ k\left(1-f_4\right)(f_2^{″}-\sigma RT{f}_5^{″})+\alpha f_4^{\prime }\left(1+f_1^{\prime }\right)+(1-f_4)q_V=0,\hfill \\ D{f}_5^{″}+kS{\left(f_5(f_2^{\prime }-\sigma RT{f}_5^{\prime })\right)}^{\prime }+\alpha {\left(f_4{f}_5\right)}^{\prime }\left(1+f_1^{\prime }\right)+2\alpha f_4{f}_5{f}_1^{″}+q_S=0,\hfill \end{array}$$

and

$$\begin{array}{c}k(f_2^{″}-\sigma RT{f}_5^{″})+q_V-2{\alpha }_1=0,\hfill \\ {\lambda }^{\ast }{f}_1^{″}-(f_2^{\prime }-\sigma RT{f}_5^{\prime })-{\alpha }_1^{2}x{f}_3=0,\hfill \\ k({\rho }_F^0-f_3)(f_2^{″}-\sigma RT{f}_5^{″})-{\alpha }_{1}x{f}_3^{\prime }+({\rho }_F^0-f_3)q_V=0,\hfill \\ k\left(1-f_4\right)(f_2^{″}-\sigma RT{f}_5^{″})-{\alpha }_{1}x{f}_4^{\prime }+(1-f_4)q_V=0,\hfill \\ D{f}_5^{″}+kS{\left(f_5(f_2^{\prime }-\sigma RT{f}_5^{\prime })\right)}^{\prime }-{\alpha }_{1}x{\left(f_4{f}_5\right)}^{\prime }-2{\alpha }_1{f}_4{f}_5+q_S=0,\hfill \end{array}$$

(13)

respectively.

In the case ${\alpha }_0=0$, we may set $\alpha =1$ and derive the ansatz

$$\begin{array}{c}u=f_1\left(t\right)+\frac{{\alpha }_1}{2}{x}^2,p=f_2\left(t\right),\rho =f_3\left(t\right),\hfill \\ {\theta }_F=f_4\left(t\right),c=f_5\left(t\right),\hfill \end{array}$$

which again is not useful from the applicability point of view.

The systems of ODEs derived above are still nonlinear and their integrability is a highly nontrivial task. For example, the known handbook for exact solving ODEs [29], which seems to be the most complete book about exact solutions of ODEs, does not contains similar systems. In Section 3, we examine in detail only the ODE system system (10) with $\alpha ={\alpha }_1=0$ in order for find steady-state solutions. We were able also to construct a nontrivial solution of the ODE system (13).

It can be noted that the structures of the first, third and fourth equations in (13) are the same and it allows us to find the functions $f_3$ and $f_4$. Moreover, assuming $f_5=const=C_0$, other equations of system (13) can be integrated in a straightforward way. Finally, using ansatz (12) the following exact solution of system (1)–(5) with $\kappa =0$

$$\begin{array}{c}u=u_0+u_{1}x+{\alpha }_{1}tx+\frac{{\alpha }_1^2{\rho }_0}{{\lambda }^{\ast }}xlnx+\frac{{\alpha }_1^2{\rho }_F^0}{6{\lambda }^{\ast }}{x}^3+\frac{1}{{\lambda }^{\ast }\gamma }\left(p_2{e}^{\gamma x}-p_1{e}^{-\gamma x}\right),\hfill \\ p=p_1{e}^{-\gamma x}+p_2{e}^{\gamma x}+p_0,\hfill \\ \rho ={\rho }_F^0+\frac{{\rho }_0}{x^2},{\theta }_F=1+\frac{{\theta }_0}{x^2},c=C_0,\hfill \end{array}$$

(14)

has been derived. Here $u_0,u_1,p_1,p_2,{\rho }_0$ and ${\theta }_0$ are arbitrary constants, other parameters are

$$\begin{array}{c}{p}_0=\frac{a{l}_P{P}_B-q_0-2{\alpha }_1+a{l}_P{\sigma }_{C}RT(C_0-C_B)}{a{l}_P+q_1},\hfill \\ C_0=\frac{a{l}_P{C}_B{F}_1{S}_C}{q_1(S-1)+a{l}_P(S-F_2{S}_C)},\gamma =\sqrt{\frac{a{l}_P+q_1}{k}},\hfill \\ {\alpha }_1=\frac{1}{2(1-S)}\left(\frac{a{p}_D{C}_B}{C_0}-q_1{P}_B-a{p}_D+q_1{P}_{B}S-q_0(1-S)+q_1{\sigma }_{C}RT(1-S)(C_B-C_0)\right).\hfill \end{array}$$

Obviously, the exact solution (14) reduces to a steady-state solution provided ${\alpha }_1=0$.

It should be pointed out that functions $q_L,q_V$ and $q_S$ of system (1)–(5) can be assumed as arbitrary functions that satisfy some basic physical laws. In such a case, the problem of Lie symmetry classification (Lie group classification) is obtained, which was firstly formulated by Ovsiannikov [30]. During the last two decades, new techniques have been worked out and applied for solving this problem (see [27,31,32] and papers cited therein). We are going to treat this problem in a forthcoming paper together with exact solutions derived via Lie symmetries.

3. Steady-State Solutions of the Model for Perfused Tissue Undergoing Deformation

At the steady-state fluid and solute fluxes have to stabilize. There is also no change of the effective stress tensor. In this case, we arrive at a system consisting of three ordinary differential equations only because the 3rd and 4th equations of (1)–(5) coincide with the first equation. So, we obtain the following equations describing steady-state tissue deformation, and the corresponding changes of interstitial pressure and concentration distribution within the tissue

$$\begin{array}{c}k(P^{″}-\sigma RT{C}^{″})+a{l}_P(P_B-P-{\sigma }_{C}RT(C_B-C))-q_L=0,\hfill \\ \frac{d}{dx}\left(-(P-\sigma RTC)+{\lambda }^{\ast }{U}^{\prime }+\kappa {\left(U^{\prime }\right)}^2\right)=0,\hfill \\ D{C}^{″}+kS\frac{d}{dx}\left((P^{\prime }-\sigma RT{C}^{\prime })C\right)+a{p}_D(C_B-C)\hfill \\ +S_C{C}_{M}a{l}_P\left(P_B-P-{\sigma }_{C}RT(C_B-C)\right)-C{q}_L=0\hfill \end{array}$$

(15)

(hereafter the upper sign ${}^{\prime }$ means the derivation $\frac{d}{dx}$). In what follows we assume

$$q_L\left(x\right)=q_0+q_{1}P\left(x\right),C_M\left(x\right)=F_1{C}_B+F_{2}C\left(x\right),F_1+F_2=1.$$

Notably, the above restrictions are widely used (see, e.g., [33]).

Let us introduce the effective pressure

$$h\left(x\right)=P\left(x\right)-\sigma RTC\left(x\right)$$

(16)

(i.e., h is a new unknown function) in order to simplify the ODE system for finding steady-states. So, we arrive at the ODE system

$$\begin{array}{c}k{h}^{″}-\left(a{l}_P+q_1\right)h-RT\left(a{l}_P\left(\sigma -{\sigma }_C\right)+q_1\sigma \right)C+a{l}_P(P_B-C_{B}RT{\sigma }_C)-q_0=0,\hfill \\ \frac{d}{dx}\left(-h+{\lambda }^{\ast }{U}^{\prime }+\kappa {\left(U^{\prime }\right)}^2\right)=0,\hfill \\ D{C}^{″}+kS\frac{d}{dx}\left(h^{\prime }C\right)-RT\left(a{l}_P{F}_2{S}_C\left(\sigma -{\sigma }_C\right)+q_1\sigma \right)C^2\hfill \\ -\left(a{p}_D+q_0+a{l}_P{C}_B{S}_{C}RT\left(F_1\left(\sigma -{\sigma }_C\right)+F_2{\sigma }_C\right)+a{l}_P{F}_2{S}_C\left(h-P_B\right)+q_{1}h\right)C\hfill \\ +a{p}_D{C}_B+a{l}_P{C}_B{F}_1{S}_C\left(P_B-C_{B}RT{\sigma }_C-h\right)=0.\hfill \end{array}$$

(17)

In contrast to the corresponding system neglecting the internal sources/sinks (see [21]), construction of the general solution of the nonlinear three-component system (17) is a highly nontrivial task. We were unable to construct its general solution. Here we restrict ourselves to search some particular exact solutions. It can be noted that there are two essentially different cases

$$\mathbf{1}.a{l}_P\left(\sigma -{\sigma }_C\right)+q_1\sigma \ne 0$$

and

$$\mathbf{2}.a{l}_P\left(\sigma -{\sigma }_C\right)+q_1\sigma =0,$$

(18)

because the latter leads to the autonomous equation for the function h.

3.1. Analysis of Case 1

In this case, one can express the function C via the function h and its second-order derivative using the first equation from (17). Substituting the function obtained into the last equation of (17), one arrives at a nonlinear fourth-order ODE with a very cumbersome structure. In order to avoid solving the latter, we assume that the concentration C is a positive constant to be determined therefore the first ODE of system (17) can be easily integrated and its general solution reads as

$$h=h_1{e}^{-\gamma x}+h_2{e}^{\gamma x}+h_3,h_1^2+h_2^2\ne 0,\gamma \in \mathbb{R}\backslash \left\{0\right\},$$

(19)

where $h_1$ and $h_2$ are arbitrary constants, $\gamma =\sqrt{\frac{a{l}_P+q_1}{k}}$, while $h_3$ is determined by the constant C as follows

$$h_3=-RT(\sigma -{\sigma }_C)C+(P_B-C_{B}RT{\sigma }_C)-\frac{q_0}{a{l}_p}.$$

Substituting (19) into the last ODE of system (17) with $q_1=0$, the concentration was specified as $C\left(x\right)=\frac{F_1{S}_C}{S-F_2{S}_C}{C}_B$ under the coefficient restrictions which are listed below. Finally one needs to integrate the second ODE with respect to the function $U\left(x\right)$. Thus, the following exact solution of system (17) with $q_1=0$ is derived

$$\begin{array}{c}P\left(x\right)=h_1{e}^{-\gamma x}+h_2{e}^{\gamma x}+h_0,\hfill \\ U\left(x\right)=\left\{\begin{array}{c}{u}_0-\frac{{\lambda }^{\ast }}{2\kappa }x\pm \frac{1}{\kappa }{\int }_0^x\sqrt{u_1+\kappa \left(h_1{e}^{-\gamma y}+h_2{e}^{\gamma y}\right)}dy,\mathrm{if}\kappa \ne 0,\hfill \\ u_0+u_{1}x-\frac{h_1}{\gamma {\lambda }^{\ast }}{e}^{-\gamma x}+\frac{h_2}{\gamma {\lambda }^{\ast }}{e}^{\gamma x},\mathrm{if}\kappa =0,\hfill \end{array}\right.\hfill \\ C\left(x\right)=\frac{F_1{S}_C}{S-F_2{S}_C}{C}_B,\hfill \end{array}$$

(20)

where $u_0,u_1,h_1$ and $h_2$ are arbitrary constants, while

$$\gamma =\sqrt{\frac{a{l}_P}{k}},q_0=\frac{a{p}_D(S-S_C)}{F_1(1-S)S_C},h_0=P_B-\frac{q_0}{a{l}_p}-\frac{C_{B}RT{\sigma }_C}{S-F_2{S}_C}(S-S_C).$$

Obviously the concentration $C\left(x\right)>0$ provided $S-F_2{S}_C>0$ and the latter always takes place because typically $S>S_C$. Moreover, the restriction $S>S_C$ is needed in order to have $q_0>0$.

In the case $q_1\ne 0$, the function (19) can be derived as well, however, the coefficient restrictions are too cumbersome. Therefore we restricted ourselves to the special case $S=S_C$ (i.e., $\sigma ={\sigma }_C$). Having this restriction, the exact solution of system (17) was derived in the form

$$\begin{array}{c}P\left(x\right)=h_1{e}^{-\gamma x}+h_2{e}^{\gamma x}+h_0,\hfill \\ U\left(x\right)=\left\{\begin{array}{c}{u}_0-\frac{{\lambda }^{\ast }}{2\kappa }x\pm \frac{1}{\kappa }{\int }_0^x\sqrt{u_1+\kappa \left(h_1{e}^{-\gamma y}+h_2{e}^{\gamma y}\right)}dy,\mathrm{if}\kappa \ne 0,\hfill \\ u_0+u_{1}x-\frac{h_1}{\gamma {\lambda }^{\ast }}{e}^{-\gamma x}+\frac{h_2}{\gamma {\lambda }^{\ast }}{e}^{\gamma x},\mathrm{if}\kappa =0,\hfill \end{array}\right.\hfill \\ C\left(x\right)=\frac{a{l}_p{F}_{1}S}{a{l}_p{F}_{1}S-q_1(1-S)}{C}_B,\hfill \end{array}$$

where $u_0,u_1,h_1$ and $h_2$ are arbitrary constants, while

$$\gamma =\sqrt{\frac{a{l}_P+q_1}{k}},q_0=-q_1\left(P_B+\frac{a{p}_D}{a{l}_p{F}_{1}S}+\frac{C_{B}RT\sigma q_1(1-S)}{a{l}_p{F}_{1}S-q_1(1-S)}\right),$$

$$h_0=\frac{a{l}_p}{a{l}_p+q_1}\left[P_B-\frac{q_0}{a{l}_p}+\frac{C_{B}RT\sigma q_1(1-S)}{a{l}_p{F}_{1}S-q_1(1-S)}\right].$$

Obviously $C\left(x\right)>0$ provided $q_1<\frac{a{l}_p{F}_{1}S}{1-S}$. The latter leads to the inequality $q_0<0$, which is incompatible with our assumption on free parameters constant of $q_L$.

3.2. Analysis of Case 2

In this case, the parameter $q_1$ is specified as follows $q_1=a{l}_P\frac{{\sigma }_C-\sigma }{\sigma }$. So, we again need $S>S_C$ or at least $S=S_C$.

Firstly we examine a special case when $h\left(x\right)=h_0.$ Obviously, the constant $h_0$ can be easily identified from the first equation of system (17):

$$h_0=\frac{\sigma }{{\sigma }_C}\left(P_B-\frac{q_0}{a{l}_P}-C_{B}RT{\sigma }_C\right).$$

(21)

The general solution of the third equation of system (17) takes the form

$$x=\pm \int {\left({\beta }_3{C}^3+{\beta }_2{C}^2+{\beta }_{1}C+{\beta }_0\right)}^{-\frac{1}{2}}dC+x_0,$$

where

$${\beta }_1=-\frac{2C_B}{D{\sigma }_C}\left[F_1{S}_C\left(q_0\sigma +a{l}_P\left(\sigma -{\sigma }_C\right)\left(C_{B}RT{\sigma }_C-P_B\right)\right)+a{p}_D{\sigma }_C\right],$$

$${\beta }_2=\frac{a{p}_D}{D}+\frac{(1-F_2{S}_C)}{D{\sigma }_C}\left(q_0\sigma +a{l}_P{P}_B\left({\sigma }_C-\sigma \right)\right)+\frac{a{l}_P{C}_{B}RT}{D}\left(\sigma -{\sigma }_C\right)\left(1+(F_1-F_2)S_C\right),$$

$${\beta }_3=\frac{2a{l}_{P}RT}{3D}(1-F_2{S}_C)\left({\sigma }_C-\sigma \right),$$

while ${\beta }_0$ and $x_0$ are arbitrary constants (the latter can be set zero without losing a generality). Integrating the second equation of system (17) and using (16), we arrive at the steady-state solution

$$\begin{array}{c}x=\pm \int {\left({\beta }_3{C}^3+{\beta }_2{C}^2+{\beta }_{1}C+{\beta }_0\right)}^{-\frac{1}{2}}dC,\hfill \\ P\left(x\right)=\sigma \left(\frac{P_B}{{\sigma }_C}-\frac{q_0}{a{l}_P{\sigma }_C}-C_{B}RT+RTC\left(x\right)\right),\hfill \\ U\left(x\right)=u_0+u_{1}x,\hfill \end{array}$$

(22)

of the initial system (15) (here ${\beta }_0,u_0$ and $u_1$ are arbitrary constants).

Notable, in this case the function $U\left(x\right)$ is linear independently on the value of $\kappa$.

Remark 1.

The integral in (22) leads to an elliptic function. Setting ${\beta }_0=0$ (the integral with ${\beta }_0\ne 0$ can be easily reduced to that with ${\beta }_0=0$), we obtain

$$x=\pm \frac{2}{{({{\beta }_2}^2-4{\beta }_1{\beta }_3)}^{1/4}}\mathit{EllipticF}\left(\sqrt{\frac{2{\beta }_{3}C}{{\beta }_2+\sqrt{{{\beta }_2}^2-4{\beta }_1{\beta }_3}}+1},\frac{1}{\sqrt{2}}\sqrt{\frac{{\beta }_2}{\sqrt{{{\beta }_2}^2-4{\beta }_1{\beta }_3}}+1}\right).$$

The case ${{\beta }_2}^2=4{\beta }_1{\beta }_3$ is special and the functions arctan or arctanh are obtained depending on sign of ${\beta }_2$.

Now we analyze system (17) without any restrictions on $h\left(x\right)$. Since the condition (18) takes place the general solution of the first equation in (17) can be easily derived

$$h\left(x\right)=h_1{e}^{-\gamma x}+h_2{e}^{\gamma x}+h_0,$$

(23)

where $h_1$ and $h_2$ are arbitrary constants ($h_1^2+h_2^2\ne 0$ otherwise we arrive at formulae (22)), while $\gamma =\sqrt{\frac{a{l}_P{\sigma }_C}{k\sigma }}$ and $h_0$ is defined above. Since the function $h\left(x\right)$ has the same structure as in (19), the formula for the function $U\left(x\right)$ is the same as in (20).

Now we turn to the function $C\left(x\right)$. Substituting (23) into the third equation of (17), one arrives at the ODE

$$\begin{array}{c}D{C}^{″}+kS\frac{d}{dx}\left(h^{\prime }C\right)+a{l}_{P}RT(1-F_2{S}_C)\left(\sigma -{\sigma }_C\right)C^2\hfill \\ +\left(H_2(h-h_0)+H_1\right)C-a{l}_P{C}_B{F}_1{S}_C(h-h_0)+H_0=0,\hfill \end{array}$$

(24)

where

$$H_2=a{l}_P\left(1-F_2{S}_C-\frac{{\sigma }_C}{\sigma }\right),$$

$$H_1=h_0{H}_2+a{l}_P{P}_B{F}_2{S}_C-a{p}_D-q_0-a{l}_P{C}_{B}RT{S}_C\left(F_1\left(\sigma -{\sigma }_C\right)+F_2{\sigma }_C\right)$$

and

$$H_0=a{p}_D{C}_B+a{l}_P{C}_B{F}_1{S}_C(P_B-C_{B}RT{\sigma }_C-h_0).$$

Subcase 1. In the case $\sigma -{\sigma }_C\ne 0$ Equation (24) is a nonlinear ODE with nonconstant coefficients, which is nonintegrable. However, we were able find some particular solutions of this equation using ad hoc ansatz

$$C\left(x\right)={\mu }_0+{\mu }_1{e}^{-\nu x}+{\mu }_2{e}^{\nu x}$$

with the to-be-determined parameters ${\mu }_0,{\mu }_1,{\mu }_2$ and $\nu$.

Substituting this ansatz into Equation (24), it was proved that an exact solution exists if $\nu =\frac{\gamma }{2}$ or $\nu =\gamma$, $\gamma =\sqrt{\frac{a{l}_P{\sigma }_C}{k\sigma }}$ and some additional restrictions take place.

In the case

$$C\left(x\right)={\mu }_1{e}^{-\frac{\gamma }{2}x}+{\mu }_2{e}^{\frac{\gamma }{2}x}+{\mu }_0$$

(25)

the restrictions

$$\begin{array}{c}{F}_2=\frac{2\sigma +(3S-2){\sigma }_C}{2\sigma S_C},\hfill \\ h_1=\frac{RT{\sigma }_C(2-3S)\left(\sigma -{\sigma }_C\right)}{2C_B{F}_1{S}_C\sigma +S{\sigma }_C{\mu }_0}{\mu }_1^2,h_2=\frac{RT{\sigma }_C(2-3S)\left(\sigma -{\sigma }_C\right)}{2C_B{F}_1{S}_C\sigma +S{\sigma }_C{\mu }_0}{\mu }_2^2\hfill \end{array}$$

hold, while the coefficients ${\mu }_0,{\mu }_1$ and ${\mu }_2$ are defined by the equations

$$\begin{array}{c}2(H_0+H_1{\mu }_0)\sigma +a{l}_{P}RT{\sigma }_C\left(2-3S\right)\left(\sigma -{\sigma }_C\right)\left({\mu }_0^2+2{\mu }_1{\mu }_2\right)=0,\hfill \\ H_1\sigma {\mu }_1+\frac{Da{l}_p{\sigma }_C}{4k}{\mu }_1-a{l}_{p}S{\sigma }_C{h}_1{\mu }_2+a{l}_{p}RT{\sigma }_C(2-3S)\left(\sigma -{\sigma }_C\right){\mu }_0{\mu }_1=0,\hfill \\ H_1\sigma {\mu }_2+\frac{Da{l}_p{\sigma }_C}{4k}{\mu }_2-a{l}_{p}S{\sigma }_C{h}_2{\mu }_1+a{l}_{p}RT{\sigma }_C(2-3S)\left(\sigma -{\sigma }_C\right){\mu }_0{\mu }_2=0.\hfill \end{array}$$

The function

$$C\left(x\right)={\mu }_1{e}^{-\gamma x}+{\mu }_2{e}^{\gamma x}+{\mu }_0$$

(26)

is an exact solution of Equation (24) if essentially different restrictions take place depending on values of ${\mu }_1$ and ${\mu }_2$.

If ${\mu }_1{\mu }_2\ne 0$ then

$$h_1=\frac{RT\sigma (1-F_2{S}_C)\left(\sigma -{\sigma }_C\right)}{(1-2S){\sigma }_C-(1-F_2{S}_C)\sigma }{\mu }_1,h_2=\frac{RT\sigma (1-F_2{S}_C)\left(\sigma -{\sigma }_C\right)}{(1-2S){\sigma }_C-(1-F_2{S}_C)\sigma }{\mu }_2,$$

and ${\mu }_0,{\mu }_1$ and ${\mu }_2$ should satisfy the equations

$$\begin{array}{c}{H}_0+H_1{\mu }_0+(h_2{\mu }_1+h_1{\mu }_2)H_2+a{l}_{P}RT\left(1-F_2{S}_C\right)\left(\sigma -{\sigma }_C\right)\left({\mu }_0^2+2{\mu }_1{\mu }_2\right)=0,\hfill \\ H_1{\mu }_1+H_2{h}_1{\mu }_0+2a{l}_{P}RT\left(1-F_2{S}_C\right)\left(\sigma -{\sigma }_C\right){\mu }_0{\mu }_1\hfill \\ -a{l}_P{h}_1{C}_B{F}_1{S}_C+\frac{a{l}_P{\sigma }_C}{k\sigma }\left(kS{h}_1{\mu }_0+D{\mu }_1\right)=0.\hfill \end{array}$$

(27)

If ${\mu }_1{\mu }_2=0$ then we may set ${\mu }_1\ne 0,{\mu }_2=0$ without losing a generality. As a result

$$h_1=\frac{RT\sigma (1-F_2{S}_C)\left(\sigma -{\sigma }_C\right)}{(1-2S){\sigma }_C-(1-F_2{S}_C)\sigma }{\mu }_1,$$

while ${\mu }_0,{\mu }_1$, and $h_2$ should satisfy the equations

$$\begin{array}{c}{h}_2\left(a{l}_P{C}_B{F}_{1}S\sigma -\left(H_2\sigma +a{l}_{P}S{\sigma }_C\right){\mu }_0\right)=0,\hfill \\ H_0+H_1{\mu }_0+H_2{h}_2{\mu }_1+a{l}_{P}RT\left(1-F_2{S}_C\right)\left(\sigma -{\sigma }_C\right){\mu }_0^2=0,\hfill \\ H_1{\mu }_1+H_2{h}_1{\mu }_0+2a{l}_{P}RT\left(1-F_2{S}_C\right)\left(\sigma -{\sigma }_C\right){\mu }_0{\mu }_1\hfill \\ -a{l}_P{h}_1{C}_B{F}_1{S}_C+\frac{a{l}_P{\sigma }_C}{k\sigma }\left(kS{h}_1{\mu }_0+D{\mu }_1\right)=0.\hfill \end{array}$$

(28)

Obviously the appropriate parameters ${\mu }_i$ and $h_2$ can be found by direct calculations because (27) and (28) are systems of algebraic equations.

Remark 2.

The 2nd and 3rd equations in (28) are nothing else but (27) with ${\mu }_2=0$.

Finally, using (16), we obtain the hydrostatic pressure

$$P\left(x\right)=h_1{e}^{-\gamma x}+h_2{e}^{\gamma x}+\sigma \left(\frac{P_B}{{\sigma }_C}-\frac{q_0}{a{l}_P{\sigma }_C}-C_{B}RT+RTC\left(x\right)\right),$$

where $C\left(x\right)$ is specified in (25) and (26).

Now we present a simple example setting the parameter ${\mu }_0=0$. In this special case, system (28) is transformed into

$$\begin{array}{c}{F}_1{h}_2=0,H_0+H_2{h}_2{\mu }_1=0,H_1{\mu }_1-a{l}_P{C}_B{F}_1{S}_C{h}_1+\frac{a{l}_{P}D{\sigma }_C}{k\sigma }{\mu }_1=0.\hfill \end{array}$$

So, two different subcases spring up: $h_2=0$ and $h_2\ne 0$. As a result, we arrive at the restrictions of the form either

$$\begin{array}{c}{\mu }_0={\mu }_2=h_2=0,h_1=\frac{RT\sigma (1-F_2{S}_C)\left(\sigma -{\sigma }_C\right)}{(1-2S){\sigma }_C-(1-F_2{S}_C)\sigma }{\mu }_1,\hfill \\ q_0=\frac{a{l}_P}{\sigma }\left({\sigma }_C-\sigma \right)\left(C_{B}RT{\sigma }_C-P_B\right)-\frac{a{p}_D{\sigma }_C}{F_1\sigma S_C},\hfill \\ D=\frac{k{C}_B{F}_{1}RT{S}_C\sigma \left(\sigma -{\sigma }_C\right)(1-2S)}{(1-2S){\sigma }_C-(1-F_2{S}_C)\sigma }-\frac{ka{p}_D\sigma }{a{l}_P{F}_1{S}_C},\hfill \end{array}$$

or

$$\begin{array}{c}{\mu }_0={\mu }_2=F_1=0,h_1=\frac{RT\sigma (\sigma -{\sigma }_C){\mu }_1}{1-2S-\sigma },h_2=\frac{a{p}_D{C}_B\sigma }{a{l}_{P}S{\sigma }_C{\mu }_1},\hfill \\ q_0=a{l}_P\left(P_B-C_{B}RT{\sigma }_C+\frac{D{\sigma }_C}{k{\sigma }^2}+\frac{C_{B}RT{\sigma }_C^2-P_B{\sigma }_C}{\sigma }\right)-\frac{a{p}_D}{\sigma }.\hfill \end{array}$$

The parameter ${\mu }_1$ is an arbitrary nonzero constant. Thus, the functions $P\left(x\right)$ and $C\left(x\right)$ can be explicitly written down.

Subcase 2. In the case ${\sigma }_C=\sigma$ (then automatically $S_C=S$), Equation (24) is the linear ODE with nonconstant coefficients. Note that this restriction leads to $q_1=0$ (see (18)). The general solution of Equation (24) with arbitrary coefficients $h_1$ and $h_2$ can be expressed in the terms of hypegeometric functions (see integration of a similar equation in [34]), which are inconvenient for practical applications. However, there are several cases when the general solution can be constructed in a simpler form provided some coefficient restrictions take place.

Let us assume additionally that $F_2=0$, i.e., $C_M=C_B$. In this case, taking into account (21) we obtain

$$H_0=(a{p}_D+S{q}_0)C_B,H_1=-(a{p}_D+q_0),H_2=0.$$

Thus, Equation (24) takes form

$$\begin{array}{c}D{C}^{″}+k\gamma S(-h_1{e}^{-\gamma x}+h_2{e}^{\gamma x})C^{\prime }+\left(a{l}_{P}S(h_1{e}^{-\gamma x}+h_2{e}^{\gamma x})-(a{p}_D+q_0)\right)C\hfill \\ =a{l}_{P}S(h_1{e}^{-\gamma x}+h_2{e}^{\gamma x})C_B-(a{p}_D+S{q}_0)C_B.\hfill \end{array}$$

(29)

Interestingly the above ODE possesses the simple particular solutions $C\left(x\right)=C_B$ provided $S=1$.

To avoid cumbersome formulae, we present the general solution of ODE (29) for $h_1=0,h_2\ne 0$ and under the additional restriction $D=\frac{a{p}_D+q_0}{{\gamma }^2}$. In this case, we have

$$\begin{array}{c}C\left(x\right)=C_1{e}^{-\gamma x}+C_2\left(A_1{e}^{\gamma x}+1\right)exp\left(-A_1{e}^{\gamma x}-\gamma x\right)\hfill \\ +\frac{exp\left(-A_1{e}^{\gamma x}-\gamma x\right)}{A_1}[(C_B-A_2)(A_1{e}^{\gamma x}+1)Ei(1,-A_1{e}^{\gamma x})\hfill \\ +\left(A_1{C}_B{e}^{\gamma x}+\gamma C_{B}x-\gamma A_{2}x-A_2\right)exp\left(A_1{e}^{\gamma x}\right)],\hfill \end{array}$$

where $C_1$ and $C_2$ are arbitrary constants while $\gamma =\sqrt{\frac{a{l}_P}{k}}$, $A_1=\frac{a{l}_P{h}_{2}S}{a{p}_D+q_0},A_2=\frac{a{p}_D{C}_B+q_{0}S{C}_B}{a{p}_D+q_0}$ and $Ei(1,-e^y)=\int exp\left(e^y\right)dy$ is the exponential integral.

The condition ${\sigma }_C=\sigma$ was used in our studies [34,35] devoted to the peritoneal dialysis transport, in which we assumed that Staverman reflection coefficient is the same for the tissue and blood capillary wall. Introduction of the tissue reflection coefficient $\sigma =1-S$ allows for the application of so called extended Darcy’s formula. In this approach, volumetric flux across the tissue is described as dependent on the hydrostatic pressure gradient decreased by the effective osmotic gradient (colloid or crystalloid). Therefore, the parameter $\sigma$ reflects some retardation factor, exerted on the fluid and solute transport across the tissue, related to the tissue composition. Thus, relation between S and $S_c$ might vary between tissue types and different conditions (fibrosis or inflammatory processes) reflecting retardation properties of this two barriers: blood capillary wall and tissue. In the case of modeling of small solutes transport across the abdominal wall muscle during peritoneal dialysis, the inequality $S>S_c$ is typically assumed [7] because $\sigma$ is rather small. However, much higher values of retardation factor ($\sigma =0.5$) were assumed in the case of IgG transport during intraperitoneal therapy [9].

4. Examples of Applications to Real-World Processes

The families of exact solutions constructed in the previous section involve some free parameters (arbitrary constants). Now we consider the basic Equation (17) with some boundary conditions, which reflect the relevant real-world situations.

In many experiments on poroelastic media such as tissue or cartilage properties, the external force is applied to obtain given shrinkage or stretching. However, such experiments are performed for non perfused material, and were investigated previously [21]. Due to the lack of similar experiments looking at the impact of the blood perfusion on the deformation properties of the tissue, in the following examples we present purely numerical predictions/simulations. In the following examples we are looking for a hydrostatic pressure, solute concentration, and tissue deformation distribution of a homogeneous material (e.g., the perfused tissue) that due to some external forces, applied to obtain assumed overall shrinkage or stretching, evolves to a new equilibrium. We look at the solutions for a fixed deformation of tissue with various elastic properties and under various boundary conditions.

Let us assume that initially $U=0$ for all x and the initial tissue width $L_0=1.0$ cm after loading external force has reached width $L_S=1.1$ cm of the deformed material at a new equilibrium state. We also assume that tissue left boundary is fixed at $x=0$ cm, whereas remaining part of the tissue is moving upon the external force. For both examples hydraulic conductivity of the tissue was assumed equal to $k=5.15\times {10}^{-5}$${\mathrm{cm}}^2/\min /\mathrm{mmHg}$ whereas hydraulic conductance for the blood capillary wall was taken as $a{l}_P=7.3\times {10}^{-5}$ 1/min/mmHg [35].

Example 1. Let us consider example that corresponds to the experiments with external pressure of $P_{ex}$ and deformation of $U=L_S-L_0$. Such model can be described by the basic Equations (15) with boundary conditions

$$P\left(0\right)=P_0,P\left(L_S\right)=P_{ex},U\left(0\right)=0,U\left(L_S\right)=L_S-L_0,C\left(0\right)=C_0,C\left(L_S\right)=C_{ex},$$

(30)

where all parameters arising in RHS are assumed some known positive constants. Mathematically it means that we consider the Dirichlet boundary conditions. Let us specify the parameters in the exact solution (20), which guarantee that formulae (20) give the exact solution of BVP (15) and (30).

Obviously, the parameters for the concentration should satisfy the restriction

$$C_0=C_{ex}=\frac{F_1{S}_C}{S-F_2{S}_C}{C}_B,$$

(31)

otherwise the solution in question does not exists. Direct calculations show that solution (20) satisfies conditions (30) if the constants $u_0,u_1,h_1$ and $h_2$ have the forms

$$\begin{array}{c}{h}_1=-h_0+\frac{h_0}{1+exp\left(\gamma L_S\right)}-\frac{P_{0}exp\left(2\gamma L_S\right)-P_{ex}exp\left(\gamma L_S\right)}{1-exp\left(2\gamma L_S\right)},h_2=-\frac{h_0}{1+exp\left(\gamma L_S\right)}+\frac{P_0-P_{ex}exp\left(\gamma L_S\right)}{1-exp\left(2\gamma L_S\right)},\hfill \\ u_0=\left\{\begin{array}{c}0,\mathrm{if}\kappa \ne 0,\hfill \\ \frac{h_1-h_2}{\gamma {\lambda }^{\ast }},\mathrm{if}\kappa =0,\hfill \end{array}\right.\hfill \\ \left\{\begin{array}{c}\pm {\int }_0^{L_S}\sqrt{u_1+\kappa \left(h_1{e}^{-\gamma y}+h_2{e}^{\gamma y}\right)}dy=\kappa \left(L_S-L_0\right)+\frac{{\lambda }^{\ast }{L}_S}{2},\mathrm{if}\kappa \ne 0,\hfill \\ u_1=\frac{L_S-L_0}{L_S}+\frac{h_2-h_1+h_{1}exp\left(-\gamma L_S\right)-h_{2}exp\left(\gamma L_S\right)}{\gamma {\lambda }^{\ast }{L}_S},\mathrm{if}\kappa =0.\hfill \end{array}\right.\hfill \end{array}$$

Note that steady-state pressure profile P does not depend on the Terzaghi stress tensor parameters $\kappa$ and ${\lambda }^{\ast }$ but on the boundary values of hydrostatic pressures as presented in Figure 1. Interstitial pressure changes from the value $P_0$ at the left boundary to the value $P_{ex}$ at the right boundary. Moreover, as it comes from (31) solute profile across tissue remains constant depending only on the solute concentration in blood $C_B$, sieving properties of blood capillary wall and tissue $S_C$ and S, respectively, and weighting factor $F_1$, see Figure 1. In all simulations, the parameters for solute with low molecular weight were considered, such as glucose assuming $C_B=6$ mmol/L, $a{p}_D=3.4\times {10}^{-2}$ 1/min, $S_C=1-{\sigma }_C=0.450$, $S=0.451$, and weighting factor $F_1=F_2=0.5$ [35].

Figure 2 shows corresponding displacement profiles as a function of distance from $x=0$ for a given tissue deformation. On the contrary to the pressure and concentration profiles, local tissue displacement depends on the Terzaghi stress tensor parameters. In the Figure 2 (left pane), difference in the deformation profile between connective (${\lambda }^{\ast }=100$, [6]) and tumor tissue (${\lambda }^{\ast }=700$) was presented whereas impact of nonlinear part of stress tensor on the local tissue displacement is presented in Figure 2 (right panel).

Moreover, even assuming the same value of tissue stretching, i.e., having $L_S=1.1$ cm but applying various external hydrostatic pressures $P_{ex}=0,10,15,20,30$ mmHg we are getting different interstitial pressure and deformation profiles within the tissue as presented in Figure 3.

Example 2. Let us consider the basic Equations (15) with boundary conditions

$$\frac{dP}{dx}{|}_{x=0}=0,P\left(L_S\right)=P_{ex},U\left(0\right)=0,U\left(L_S\right)=L_S-L_0,\frac{dC}{dx}{|}_{x=0}=0,C\left(L_S\right)=C_{ex}.$$

(32)

Mathematically it means that we consider the zero Neumann conditions at the boundary $x=0$, i.e., no-flux conditions.

Let us specify the parameters in the exact solution (20), which guarantee that formulae (20) give the exact solution of BVP (15) and (30). Obviously, the parameters for the concentration should again satisfy the restriction (31). So, making similar calculations one shows that solution (20) satisfies conditions (32) if the constants $u_0,u_1,h_1$ and $h_2$ have the forms

$$\begin{array}{c}{h}_1=h_2=\frac{-h_0+P_{ex}}{2cosh\left(\gamma L_S\right)},u_0=0,\hfill \\ \left\{\begin{array}{c}\pm {\int }_0^{L_S}\sqrt{u_1+\kappa h_1\left(e^{-\gamma y}+e^{\gamma y}\right)}dy=\kappa \left(L_S-L_0\right)+\frac{{\lambda }^{\ast }{L}_S}{2},\mathrm{if}\kappa \ne 0,\hfill \\ u_1=\frac{L_S-L_0}{L_S}-\frac{h_0+P_{ex}}{\gamma {\lambda }^{\ast }{L}_S}tanh\left(\gamma L_S\right),\mathrm{if}\kappa =0.\hfill \end{array}\right.\hfill \end{array}$$

In the case $h_1=h_2$ and the additional restriction $u_1=2\kappa h_1$ takes place, integral in (20) can be expressed in the terms of elementary functions, hence we obtain

$$U\left(x\right)=u_0-\frac{{\lambda }^{\ast }}{2\kappa }x\pm \frac{4\sqrt{h_1}}{\gamma \sqrt{\kappa }}sinh\left(\frac{\gamma }{2}x\right).$$

Thus, the boundary condition $U\left(L_S\right)=L_S-L_0$ leads to the formula for $L_S$

$$L_S-L_0+\frac{{\lambda }^{\ast }{L}_S}{2\kappa }=\pm \frac{\sqrt{8(h+P_{ex})}}{\gamma \sqrt{\kappa }}\frac{sinh\frac{\gamma L_S}{2}}{\sqrt{cosh\left(\gamma L_S\right)}}.$$

Due to the Neumann boundary condition at the left boundary any changes of the $P_{ex}$ would lead not only to the substantial changes in the interstitial pressure profiles (similarly to the Example 1) but also to those in P in the point $x=0$, whereas the related changes in the deformation profiles would be relatively small, see Figure 4.

Example 3. Let us consider in contrast to Example 1, a small shrinkage of tissue from $L_0=1$ cm to $L_S=0.99$ cm assuming that pressures at the boundaries are the same as in Example 1. Figure 5 shows that in case of more elastic properties of the material (${\lambda }^{\ast }=100$), the impact of the boundary conditions (higher pressure closely to $x=L_S$, as in Figure 1, left profile but for $L_S=0.99$ cm instead of $L_S=1.1$ cm) on the local tissue deformation U is more pronounced. Further increase of tissue shrinkage to $L_S=0.90$ cm leads to smaller difference in the tissue displacement between considered materials with different elastic properties, Figure 5, right panel. In both cases, the relevant simulations have also showed that the impact of nonlinear term of stress tensor, $\kappa$, on the local displacement or pressure profile was negligible.

5. Discussion

Notably, our previous analysis considered nonperfused material and linear theory of poroelasticity [21]. The model presented here describe transport of fluid and solutes through the perfused tissue providing information on the changes in the hydrostatic pressure and solute concentration profiles in the tissue, and tissue deformation with a nonlinear term in the Terzaghi stress tensor [2,21,23]. The model for perfused tissue, which takes into account the local exchange between tissue and blood and lymphatic systems, is based on the system of nonlinear PDEs. Since the latter is nonintegrable, we applied the method of Lie symmetry analysis in order to identify its properties. It was established that the system possesses nontrivial properties depending on values of densities of the fluxes $q_V$, $q_L$ and the solute flux $q_S$ (Section 2). Here we used the simplest Lie symmetry that allows us to construct steady-state solutions, which are very important for the model in question. The relevant nonlinear ODE system was treated in details in order to provide closed formulae for steady-state solutions. In particular, the effect of nonlinear term in the Terzaghi stress tensor was studied (Section 3). Notably, some additional restrictions on the parameters were needed, otherwise steady-state solutions cannot be derived explicitly. Interestingly that the restriction $\sigma ={\sigma }_C$ was also used in our earlier papers [34,35], where the model has an essentially different structure.

The steady-state solutions obtained allowed us to present two realistic examples. An analysis of the impact of the parameter values and boundary conditions on the distribution of hydrostatic pressure, osmotic agent and displacement in the perfused tissue was provided (Section 4). In particular, the effect of nonlinear term in the Terzaghi stress tensor on the displacement may be evaluated (Figure 2). Thus, the solutions derived provide the insight into the role of some model parameters in the described transport and displacement processes, assuming some simplifications in the description of those physiologically and mathematically complex problems that need in general numerical approach and specification of all model parameters.

In a forthcoming paper we plan to construct time-dependent solutions of the model using the Lie symmetries identified above and to provide more complicated examples of their applications.