The Human Meniscus Behaves as a Functionally Graded Fractional Porous Medium under Confined Compression Conditions

Abstract

In this study, we observe that the poromechanical parameters in human meniscus vary spatially throughout the tissue. The response is anisotropic and the porosity is functionally graded. To draw these conclusions, we measured the anisotropic permeability and the “aggregate modulus” of the tissue, i.e., the stiffness of the material at equilibrium, after the interstitial fluid has ceased flowing. We estimated those parameters within the central portion of the meniscus in three directions (i.e., vertical, radial and circumferential) by fitting an enhanced model on stress relation confined compression tests. We noticed that a classical biphasic model was not sufficient to reproduce the observed experimental behaviour. We propose a poroelastic model based on the assumption that the fluid flow inside the human meniscus is described by a fractional porous medium equation analogous to Darcy’s law, which involves fractional operators. The fluid flux is then time-dependent for a constant applied pressure gradient (in contrast with the classical Darcy’s law, which describes a time independent fluid flux relation). We show that a fractional poroelastic model is well-suited to describe the flow within the meniscus and to identify the associated parameters (i.e., the order of the time derivative and the permeability). The results indicate that mean values of ${\lambda }_{\beta },\beta$ in the central body are ${\lambda }_{\beta }=5.5443\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{1-\beta }}$, $\beta =0.0434$, while, in the posterior and anterior regions, are ${\lambda }_{\beta }=2.851\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{1-\beta }}$, $\beta =0.0326$ and ${\lambda }_{\beta }=1.2636\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{1-\beta }}$, $\beta =0.0232$, respectively. Furthermore, numerical simulations show that the fluid flux diffusion is facilitated in the central part of the meniscus and hindered in the posterior and anterior regions.

1. Introduction

Human meniscus plays a key role in the functioning of the knee joint (Figure 1a). This tissue has a number of functions, such as: load bearing (about 45–75% of the total load on the joint), joint stability and lubrication [1,2]. Degenerative processes of the meniscus, either from injuries or ageing, affect approximately 35% of the population [3]. When required, the most common surgical procedure is currently total/partial meniscectomy (i.e., removal of the damaged tissue).

It has been noticed that, although the knee joint can still function with the total/partial absence of the meniscus, the loss of the meniscus leads to an increase in the contact stresses on the tibial/femoral cartilage proportional to the amount of meniscal tissue removed [4]. The increase in contact stresses on the cartilage is one of the main factors of Osteoarthritis (OA). Therefore, partial/total meniscal replacements are expected to help avoid articular cartilage degeneration.

Currently, the clinical outcomes of these implants are not ideal [5] due to the fact that they do not mimic the structure–property relationships of the tissue as these are still not well understood [6]. Among many biologic tissues, the meniscal tissue is composed of porous solid matrices—mainly collagen—with fluid filling the pores [7,8,9,10]. The overall mechanical behaviour of this type of tissues depends not only on the solid matrix deformation, but also on the movement of the fluid in and out of the collagen channels [7,8] during the deformation.

Meniscal tissue exhibit a non-uniform and anisotropic porosity, which is related to the graded material properties. Such graded material properties are fundamental for the correct functioning of this tissue [8]. The investigation of the material properties, coupled with the quantification of architectural parameters, such as the porosity and channel interconnectivity, will enable the design of artificial cellular structure, which can resemble the native behaviour of the tissue (see, for instance, [11]). More in depth, advanced microscopy investigations highlighted that collagen bundles form the wall of channels, which can be observed both at the macroscale and at the microscale [7,10]. Fluid is able to flow inside these channels in response to physiological loading.

Understanding and appropriately modelling the fluid flow behaviour in the different portions of the meniscal tissue, considering a range of various loading conditions, is essential to gain insight into the biomechanical function of this tissue. To date, limited information is available on region-specific and anisotropic permeability in the meniscus, which is essential to understand the fluid flow evolution and its relationship with the internal architecture.

Poromechanics experimental tests, such as confined compression tests—including both stress relaxation and creep—are currently used in order to characterize material parameters, such as the elastic modulus at equilibrium and permeability. In order to identify these two parameters, the main models used are based on the biphasic and poroelastic theories. Even though these theories developed from different roots, it can be shown that they are basically equivalent [12].

In this study, we focused our attention on the evaluation of the permeability in the meniscus and in particular on its variation in the different regions of the tissue and directions (See Figure 1b). Moreover, we showed that the biphasic model does not provide a good fit with the experimental curves; therefore, we proposed a poroelastic model in which the pore pressure diffusion equation is derived by adopting a modified version of Darcy’s law involving fractional derivatives [13].

It has been shown that, in a high porosity medium, there is a departure from Darcy’s law as the inertia (velocity-squared term), thermal dispersion, convective (development term) and boundary (no-slip condition) effects not included in the Darcy’s model may play a significant role [14]. Moreover, in the last few decades, experimental evidence of anomalous diffusive phenomena, i.e., not following a Darcian behaviour, has grown [15,16,17,18,19].

This is mainly due to the fact that the permeability, hence the rate of fluid flow, is not a constant quantity. Variations in permeability occur, for example, when the fluid flow impacts the geometry or the micro-structural features, such as the configuration of the pores. For example, experiments on water flow in building materials highlighted that the permeability changes during the flowing process as a result of the microstructural rearrangement of grains/pores.

Iaffaldano et al. [19] hinted that, during the compaction of sand, the permeability might decrease due to the fact that the fluid carries solid particles, which then close some of the pores. Essentially, the configuration of the medium, in particular the ratio between closed/open pores, changes during the process. Fluid might be trapped in the medium leading to a slower fluid flow rate. On the contrary, if during the fluid diffusion process, some of the pores open creating conductive microchannels, the permeability might increase. Therefore, fluid can be transported for a large distance in a reduced time resulting in a faster diffusion process.

Modelling the anomalous fluid diffusion process is one of the key points when dealing with the poromechanics of biological tissues. Therefore, it is fundamental to develop a theory that can incorporate the change in microstructural features (for instance, the interaction between fluid particles and open pores/channels) during the transport process. Recently, the stochastic Continuous Time Random Walk (CTRW) framework was proposed for this purpose [20]. However, a deterministic derivation, dual to the CTRW, can be represented by the introduction of a modified version of Darcy’s law involving linear fractional operators [21,22].

The goal of this paper is to present the results of experimental confined compression tests performed on samples extracted from the three portions of the human knee meniscus (posterior, central and posterior). During confined compression, fluid flows through the collagen channels with a rate depending on the permeability of the tissue itself (Figure 1c).

We aim at: (1) extracting the diagonal terms of the permeability tensor; (2) studying, analysing and modelling the observed decay of fluid flow during the test, which is not captured by the classical Darcy’s law. Hence, we propose a generalization of this law, which involves time derivatives of non integer order appropriate to model the fluid flow in the meniscal tissue; (3) proposing a small deformation fractional poroelastic model for the human meniscal tissue during confined compression. The model enables us to identify the two parameters involved in the fractional pore pressure diffusion equation, namely the permeability and the order of the (fractional) derivative.

In this study, we do not consider coupling between the flow in different directions, i.e., assume that the permeability tensor is diagonal. The structure of the paper is as follows: we introduce the rationale behind a fractional Darcy’s law, we then summarize the main equations of both biphasic and fractional poroelastic theories; we then present the confined compression poromechanics experimental tests and discuss the material parameters (fractional permeability and order of the fractional derivative) we recover through the fittings.

2. Fluid Flow in Complex Porous Media: Time-Fractional Darcy’s Law

Fluid flow in porous media is commonly modelled by Darcy’s law. This relation was observed experimentally by Darcy [23] and derived in 1986 from the Navier–Stokes equations using the homogenisation theory [24]. Darcy’s law states that the instantaneous flow rate through a homogeneously permeable porous medium of permeability k is proportional to the dynamic viscosity of the fluid and the pressure drop over a given distance.

The total discharge, $J_f$ (units $L^3/T$, where L and T indicate units of length and time, respectively) equates the product of the intrinsic permeability, k ($L^2$) with the cross-sectional area flow, A ($L^2$) and the total pressure drop $p_{\mathrm{out}}-p_{\mathrm{in}}$ ($F/L^2$), all divided by the dynamic viscosity, $\mu$ ($FT/L^2$) and the length over which the pressure drop is taking place ℓ:

$$J_f=-\frac{k·A·(p_{\mathrm{out}}-p_{\mathrm{in}})}{\ell \mu }$$

More generally, Darcy’s law states that the discharge per unit area, $j_f=J_f/A$ is proportional to the pressure gradient, and the intrinsic permeability and inversely proportional to the dynamic viscosity:

$${\mathit{j}}_{\mathit{f}}=-\frac{k}{\mu }\nabla p.$$

The derivation of Darcy’s law from Navier–Stokes equations assumes a creeping, laminar, stationary and incompressible flow of density $\rho$ and velocity ${\left(u_i\right)}_{1\le i\le 3}$. Incompressibility implies: ($\frac{D\rho u_i}{Dt}=0$), which leads to Stokes’ equation in the presence of gravity $g_i$:

$$\mu {\nabla }^2{u}_i+\rho g_i-p_{,i}=0.$$

Assuming that the viscous resisting force varies linearly with the velocity, introducing the porosity $\varphi$, and the second order permeability tensor $\mathit{k}=k_{ij}$, a simple derivation leads to:

$$-{\left(k_{ij}\right)}^{-1}\mu \varphi u_j+\rho g_i-p_{,i}=0.$$

The discharge per unit area in direction n can then be written as follows:

$$j_{fn}=-\frac{k_{ni}}{\mu }\left(p_{,i}-\rho g_i\right).$$

$${\mathit{j}}_{\mathit{f}}=-\frac{\mathit{k}}{\mu }\left(\nabla p-\rho \mathit{g}\right).$$

If the resistance to fluid flow offered by the pores varies in space, the components of the permeability tensor $k_{ij}$ vary in space. $k_{ij}$ is a symmetric (Onsager reciprocal relations), positive definite (because the flow component parallel to the pressure drop occurs in the same direction as the pressure drop) matrix. $k_{ij}$ may be isotropic, in which case, it is diagonal, and all diagonal entries are identical: $k_{ij}=k{\delta }_{ij}$. In general, the permeability tensor is anisotropic, and may also not be diagonal. In all cases, the permeability tensor can be diagonalised, as it is symmetric positive definite.

As previously underlined, when dealing with soft tissue and complex porous media, permeability changes during the deformation process; therefore, it is important to develop a theory, such as the CTRW, which incorporates the change in microstructural features during the transport process. A modified version of Darcy’s law involving linear fractional operators can be seen as the deterministic version of CTRW [21,22]:

$${\mathit{j}}_{\mathit{f}}=-\frac{\mathit{k}}{\mu }{D}_0^{\beta }\left(\nabla p\right)$$

(1)

In the following, $\frac{\mathit{k}}{\mu }$ will be indicated as ${\mathit{\lambda }}_{\mathit{\beta }}$. $D_0^{\beta }$ indicates Caputo’s fractional derivative of order $\beta$ of $\nabla p$, which is defined below:

$$(D_0^{\beta }\nabla p)\left(t\right)=\frac{1}{\mathrm{\Gamma }(n-\beta )}{\int }_a^t\frac{\nabla p^n\left(\tau \right)}{{(t-\tau )}^{\beta +1-n}}d\tau$$

(2)

The expression is valid for $n-1<\beta <n$, and $\mathrm{\Gamma }$ is the Euler’s Gamma function. In the cases considered in this paper, the lower bond of integration $a=0$. The fractional derivative method offers the possibility to model, with reduced number of parameters, all of the anomalous diffusion behaviours by changing the order of the derivative. The main drawback is that it is difficult to link the order of the derivatives with the microstructural features. Within this frame, in this paper, we identify a simple mathematical model able to describe fluid flow in the human meniscus.

3. Biphasic and Linear Fractional Poroelastic Models

The assumptions at the basis of both the poroelasticity and biphasic theories are the following:

• The solid phase is incompressible, linear elastic, subject to infinitesimal strain, homogeneous, isotropic (i.e., the material parameters of the solid phase do not depend on the orientation nor the position in the sample) and non dissipative.
• The fluid phase is incompressible, it flows slowly through the pores, it is homogeneous, isotropic and non dissipative, and there is no fluid source.
• The only dissipation comes from the frictional drag due to the relative velocities of the two phases.
• The absence of external body forces (other than those explicitly mentioned for the confined compression tests).
• The isotropy and homogeneity of the permeability tensor, which is then reduced to the scalar k representing the averaged intrinsic permeability of the sample.
• The permeability k is a constant parameter (time-independent) in the biphasic model and it is a time-dependent quantity in the fractional poroelastic model (due to the time-fractional Darcy’s law described in Equation (1)).

3.1. Biphasic Model—Consolidation Problem

For the purposes of this study, we restrict our formulation to small strain theory and consider the linear biphasic model from [25]. We assume the solid matrix is incompressible, linear elastic, isotropic, homogeneous and non dissipative, whereas the interstitial fluid is incompressible and non dissipative. The only dissipation comes from the frictional drag due to the relative velocities of the two phases. We denote by ${\lambda }_0=k/\mu$, the averaged axial permeability of the sample and $H_A$ as the aggregate modulus.

Under these assumptions and the additional assumptions coming from the confined compression tests setting, the biphasic theory leads to the following unidimensional boundary value problem for the vertical displacement $u_z^s$ of the solid phase

$$\frac{{\partial }^2{u}_z^s}{\partial z^2}=\frac{1}{H_A{\lambda }_0}\frac{\partial u_z^s}{\partial t}\mathrm{in}(0,h)\times (0,T),$$

(3)

with the boundary condition at $z=0$ where the solid skeleton is fixed

$$u_z^s(z=0,t)=0\mathrm{on}(0,T),$$

(4)

the initial condition

$$u_z^s(z,t=0)=0\mathrm{on}(0,h),$$

(5)

and finally the boundary condition at $z=h$, which will be different for the creep test and for the stress relaxation test [26]. For the creep test, the condition is

$${\frac{\partial u_z^s}{\partial z}}_{|z=h}=-\frac{P_A}{H_A}\mathrm{on}(0,T),$$

(6)

where $P_A$ is the applied compressive stress, and, for the stress relaxation test, the condition is

$$u(z=h,t)=\left\{\begin{array}{cc}-V_{0}t\hfill & \mathrm{if}0⩽t<t_0,\hfill \\ -V_0{t}_0\hfill & \mathrm{if}{t}_0⩽t⩽T,\hfill \end{array}\right.$$

(7)

where $V_0$ and $t_0$ are the input data from the stress relaxation test.

If we denote $c_n={(-1)}^{\frac{n-1}{2}}{\left(\frac{2h}{n\pi }\right)}^2$, the solution to Equations (3)–(6), in the case of the creep test, reads

$$u_z(z,t)=\frac{P_A}{H_A}\left[-z+\frac{2}{h}\sum _{n=1,3}^{\infty }{c}_{n}exp\left(-\frac{H_A{\lambda }_0{n}^2{\pi }^{2}t}{4h^2}\right)sin\left(\frac{n\pi z}{2h}\right)\right],$$

(8)

If we denote ${\rho }_n={\left(\frac{n\pi }{h}\right)}^2$, the solution to Equations (3)–(5) and (7) [26] is given by

$$u_z^s(z,t)=\left\{\begin{array}{cc}{\displaystyle -{\displaystyle \frac{V_{0}tz}{h}}-{\displaystyle \frac{2V_0}{H_A{\lambda }_{0}h}}\sum _{n=1}^{\infty }{\displaystyle \frac{{(-1)}^n}{{\rho }_n^{3/2}}}\left(1-exp(-H_A{\rho }_n{\lambda }_{0}t)\right)sin\left(\sqrt{{\rho }_n}z\right)}\hfill & \mathrm{if}0⩽t<t_0,\hfill \\ {\displaystyle -{\displaystyle \frac{V_0{t}_0}{h}}-{\displaystyle \frac{2V_0}{H_A{\lambda }_{0}h}}\sum _{n=1}^{\infty }{\displaystyle \frac{{(-1)}^n}{{\rho }_n^{3/2}}}exp(-H_A{\rho }_n{\lambda }_{0}t)\left(exp(-H_A{\rho }_n{\lambda }_{0}t)-1\right)sin\left(\sqrt{{\rho }_n}z\right)}\hfill & \mathrm{if}{t}_0⩽t⩽T.\hfill \end{array}\right.$$

(9)

The stress relaxation response ${\sigma }_t$ is given by $H_A{\frac{\partial u}{\partial z}}_{|z=h}$, and thus

$${\sigma }_t=\left\{\begin{array}{cc}{\displaystyle -{\displaystyle \frac{V_0{H}_{A}t}{h}}+{\displaystyle \frac{2V_0}{{\lambda }_{0}h}}\sum _{n=1}^{\infty }{\rho }_n^{-1}\left(1-exp(-H_A{\rho }_n{\lambda }_{0}t)\right)}\hfill & \mathrm{if}0⩽t<t_0,\hfill \\ {\displaystyle {\displaystyle \frac{2V_0}{{\lambda }_{0}h}}\sum _{n=1}^{\infty }{\rho }_n^{-1}exp(-H_A{\rho }_n{\lambda }_{0}t)\left(exp(-H_A{\rho }_n{\lambda }_{0}t)-1\right)}\hfill & \mathrm{if}{t}_0⩽t⩽T.\hfill \end{array}\right.$$

(10)

3.2. Linear Fractional Poroelastic Model: Fractional Consolidation Problem

The classical linear model of transient flow and deformation of a homogeneous fully saturated elastic porous medium depends on an appropriate coupling of the fluid pressure and solid stress. A change in applied stress produces a change in the fluid pressure or fluid mass, and a change in fluid pressure or fluid mass is responsible for a change in the volume of the porous material. The coupling term affects only the hydrostatic part of the stress tensor. The stress tensor can be written as follows:

$$\mathit{\sigma }=2G\mathit{ϵ}+\lambda \mathrm{trace}\left(\mathit{ϵ}\right)\mathbf{I}-\alpha p\mathbf{I}$$

(11)

where $\lambda =K-\frac{2}{3G}$ is the Lamé constant; $G,K$ are the shear and bulk modulus, respectively; p is the pore pressure; and $\alpha$ is the Biot coefficient. In order to solve this one, we need an additional equation, which is given by the pore pressure diffusion equation derived considering a time fractional Darcy’s law in Equation (1) [13]:

$$\frac{\partial p}{\partial t}=\frac{KB}{\alpha }{\lambda }_{\beta }{D}_0^{\beta }{\nabla }^{2}p-B\frac{\partial {\sigma }_H}{\partial t}$$

(12)

where $D_0^{\beta }$ indicates the Caputo’s time fractional derivative [27]. It is important to note that the classical (non fractional) pore pressure diffusion equation is recovered in case of $\beta =0$. In this case ${\lambda }_{\beta }=\frac{k}{\mu }$ with dimension of $\frac{{\left[L\right]}^4}{\left[F\right]\left[T\right]}$. In the case of $\beta \ne 0$ note that ${\lambda }_{\beta }$ has dimension of $\frac{{\left[L\right]}^4}{\left[F\right]{\left[T\right]}^{1-\beta }}$. The pore pressure diffusion equation can also be written in terms of strains [13]:

$$\frac{\partial p}{\partial t}=\frac{K_u-K}{{\alpha }^2}{\lambda }_{\beta }{D}_0^{\beta }{\nabla }^{2}p-\frac{K_u-K}{\alpha }\frac{\partial {ϵ}_d}{\partial t}$$

(13)

The consolidation problem is modelled through the 1D uniaxial strain poroelastic problem [13]. Equation (11) is adapted for the 1D case in which the only non-zero component of strain is ${ϵ}_{zz}$, we obtain:

$$\left(K+\frac{4G}{3}\right)\frac{\partial {ϵ}_{zz}}{\partial z}-\alpha \frac{\partial p}{\partial z}=0$$

(14)

Equation (14) is then coupled with the pore pressure diffusion equation in Equation (13) in order to obtain the following pore pressure diffusion equation:

$$\frac{\partial p}{\partial t}=\overline{\lambda }{D}_0^{\beta }\frac{{\partial }^{2}p}{\partial z^2}$$

(15)

where $\overline{\lambda }={\lambda }_{\beta }\frac{\left(4G+3K\right)\left(K_u-K\right)}{{\alpha }^2\left(4G+3K_u\right)}$. We indicate with the symbol ${\lambda }_0=\frac{\left(4G+3K\right)\left(K_u-K\right)}{{\alpha }^2\left(4G+3K_u\right)}$.

The boundary value problem is given [13]:

$$\begin{array}{c}\hfill {\frac{\partial p}{\partial z}}_{z=0}=0,\end{array}$$

(16)

$$\begin{array}{c}\hfill p(z=h,t)=0,\end{array}$$

(17)

$$\begin{array}{c}\hfill u(z=h,t)=0.\end{array}$$

(18)

A constant compressive stress in the z-direction is applied to the cylinder at $z=0$:

$${\sigma }_{zz}(0,t)=-P_A$$

(19)

where $-P_A$ is the applied compressive stress. The initial pore pressure is derived for undrained conditions, i.e.,

$$p(z,0)=P_A\frac{3(K_u-K)}{\alpha \left(4G+3K_u\right)}$$

(20)

The analytical solution in terms of pore pressure reads as follows:

$$p(z,t)=P_A\gamma \sum _{n=1,3}^{\infty }{E}_{1-\beta ,1}\left(-\frac{n^2{\pi }^2\overline{\lambda }{t}^{1-\beta }}{4h^2}\right)c_{n}cos\frac{n\pi z}{2h}$$

(21)

where:

$$\begin{array}{c}\hfill \gamma =\frac{3\left(K_u-K\right)}{\alpha \left(4G+3K_u\right)}\end{array}$$

(22)

$$\begin{array}{c}\hfill c_n={(-1)}^{\frac{n-1}{2}}{\left(\frac{2h}{n\pi }\right)}^2\end{array}$$

(23)

where $E_{1-\beta ,1}$ is the Mittag–Leffler function. In the case of $\beta =0$, the solution is identical to the classical Terzaghi’s solution in which $E_{1-\beta ,1}=exp$.

In the case of $\beta =0$, it is possible to derive the following displacement analytical solution [28]:

$$u_z(z,t)=\frac{P_A}{{\lambda }_0}\left[-z+\frac{2}{h}\sum _{n=1,3}^{\infty }{c}_{n}exp\left(-\frac{{\lambda }_0\kappa n^2{\pi }^{2}t}{4h^2}\right)sin\left(\frac{n\pi z}{2h}\right)\right],$$

(24)

where $c_n={(-1)}^{\frac{n-1}{2}}{\left(\frac{2h}{n\pi }\right)}^2$.

3.3. Correspondence of Parameters between the Biphasic and the Linear Fractional Poroelastic Model

The biphasic model essentially depends on two parameters: the aggregate modulus $H_A$ and the averaged axial permeability $\kappa$. The linear fractional poroelasticity model of Biot depends instead on three parameters: the diffusion coefficient ${\lambda }_{\beta }$, the fractional power $\beta$ and the ratio $k/\mu$. These parameters depend on the material parameters of the constituent phases or have to be determined, for example using numerical experiments curves fitting. The biphasic and linear fractional Biot models are equivalent under the above assumptions when $\beta =0$ (i.e., in the case of a non-fractional Biot’s model) and when $H_A={\lambda }_0$.

The Equations (8) and (24) show that the biphasic and Biot models are equivalent when $H_A={\lambda }_0$. Similarly, the displacement solution from the Biot’s model in the case of the stress relaxation test can be obtained by replacing $H_A$ by ${\lambda }_0$ in Equation (9). Hereinafter,

Table 1. Linear fractional Biot model.

Notation	Name	Expression	Unit
K	drained bulk modulus		Pa
G	drained shear modulus		Pa
$K_u$	undrained bulk modulus		Pa
$\alpha$	Biot’s coefficient		Unitless
B	Skempton’s coefficient	$(K_u-K)/\left(\alpha K_u\right)$	Unitless
k	averaged intrinsic permeability		${\mathrm{m}}^2$
$\mu$	fluid viscosity		$\mathrm{N}·\mathrm{s}·{\mathrm{m}}^{-2}$
$\kappa$	averaged axial permeability	$\kappa =k/\mu$	${\mathrm{m}}^4·{\mathrm{N}}^{-1}·{\mathrm{s}}^{-1}$
	non-fractional diffusion coef
${\lambda }_0$	- compressible case	$\left((4G+3K)(K_u-K)\right)/\left({\alpha }^2(4G+3K_u)\right)$	${\mathrm{m}}^3·\mathrm{s}·{\mathrm{kg}}^{-1}$
	- incompressible case	$(4G+3K)/3$
${\lambda }_{\beta }$	fractional diffusion coef.		${\mathrm{m}}^3·{\mathrm{s}}^{1+\beta }·{\mathrm{kg}}^{-1}$
$\beta$	fractional power		Unitless

and

Table 2. Biphasic model.

Notation	Name	Expression	Unit
$\gamma$	ratio solid volume/fluid volume		Unitless
$K_s$	solid phase bulk modulus		Pa
$G_s$	solid phase shear modulus		Pa
$\mathbf{k}$	intrinsic permeability		${\mathrm{m}}^2$
$\mu$	fluid viscosity		$\mathrm{N}·\mathrm{s}·{\mathrm{m}}^{-2}$
$\kappa$	averaged axial permeability	$\mu /\left({(1+\gamma )}^{2}k\right)$	$\mathrm{N}·\mathrm{s}·{\mathrm{m}}^{-4}$
$H_A$	aggregate modulus	$(4G_s+3K_s)/3$	Pa

summarize the relationships between the parameters for each models.

4. Materials and Methods

Poromechanics Tests—Confined Compression

Menisci were harvested from patients (age 65–76, mean 72, and standard deviation 4) undergoing total knee arthroplasty (ethical approval EM 249-2018 21/2017/Sper/IOR EM2, Rizzoli Orthopaedic Institute, Bologna, Italy). Samples labelled as “degraded” by gross investigation of the surgeon were discarded. Three lateral menisci and three medial menisci were collected and stored at −20 °C [29,30]. The day of the test, each meniscus was thawed in a phosphate-buffered saline (PBS) bath at room temperature for about thirty minutes [30]. Then, cylindrical samples (diameter of 3 mm and height of 3–4 mm) were extracted from the central body, anterior and posterior region along three reference directions, i.e., vertical, radial and circumferential (Figure 1b), following a dedicated procedure [31], adapted to the meniscal configuration.

A total of 18 cylindrical samples were grouped considering the the harvesting region (body, anterior and posterior) and the direction (vertical, radial and circumferential) and then tested. The testing protocol was implemented on a multi-axis mechanical tester (Mach-1, Biomomentum Inc., Laval, QC, Canada) in a confined compression configuration Figure 1c); with this specific setup, we were confident that the fluid could flow only from the base of the cylinder. Insertion in the confining chamber, thickness measurement, and removal of the meniscal sample followed a dedicated procedure [31].

Concerning confined compression, the setup involved a confining chamber with an inner diameter equal to that of the tool used during the extraction of cylindrical samples, i.e., 3 mm. The bottom of the chamber consisted of a porous-permeable platen, while the top allowed the insertion of the piston. Both these components (Biomomentum Inc., Canada) were manufactured to allow the leakage of the fluid only through the porous platen.

Moreover, an additional control was provided by visual inspection of the confining chamber transparent wall, which allowed us to check any fluid flow toward the piston (Figure 2a). We measured the amount of fluid flowing out each sample as a result of the applied compression loading by monitoring the weight of the samples throughout the test.

Testing Protocols

Relaxation tests. For the relaxation tests, we followed the procedure for permeability and aggregate modulus $H_A$ analysis, as recommended by the manufacturer of the testing machine [31,32]. In further detail, the test sequence considers a pre-compression with ramp amplitude 10% of the sample thickness (h) and ramp velocity 0.3% h/s, followed by five stress-relaxations with incremental ramp amplitude 2% h and ramp velocity 0.3% h/ s.

Creep Tests The testing protocol—specifically the loading amplitude and the detection of the fluid flow out of the sample—was designed and optimized for this specific study. In particular, three testing phases were implemented for each meniscal sample. The first phase, developed according to [31], consisted in five separated repetitions of confined compression, 75 s each, to realize 450 s of total creep. More in detail, before and after each compression, sample was removed from the confining chamber, weighed by a microbalance (Tecnopound, Ravenna, Italy) and repositioned inside the chamber with the previous vertical orientation.

The second phase involved the resting of the sample in PBS until its height returned to the pre-loading value, to thus recover the loading history. In the final phase, i.e., the third one, a single step of confined compression creep was applied to the sample for 450 s. Similarly to what we reported for the phase one, also in this case, the weight of the sample was measured before and after the compression.

It is crucial to emphasize that this third phase served to calibrate the weights measured during phase one, when the multiple removals of the sample from the confining chamber and the separated compressive ramps could have affected the creep. The stage velocity was 0.3% h, while the load target was 0.5 N, corresponding to a stress of about 0.07 MPa, which is within the range of physiological values for human menisci [33]. The decrease in weight of the sample was monitored throughout the test, and this was correlated to the amount of fluid discharged by the samples, hence, it served to reveal information regarding the rate of the fluid flow.

In particular, the weight of the sample over time $W\left(t\right)$ is related to the initial weight of the sample $W_0$ and the fluid flux:

$$W\left(t\right)=W_0-{\int }_0^t{\mathit{j}}_f·{\mathit{w}}_{s}Adt$$

(25)

with $j_f$ given by Equation (1), in which the gradient of the pressure is given by differentiating the expression of the pressure in Equation (21), $w_s$ being the specific weight and A the cross-sectional area. Substituting Equations (1) and (21) into Equation (25), we obtain the following relation:

$$W\left(t\right)=W_0-P_A\gamma \frac{2}{h}\sum _{n=1,3}^{\infty }\left[{\lambda }_{\beta }{t}^{1-\beta }{E}_{1-\beta ,2-\beta }\left(-\frac{{\pi }^2\overline{\lambda }{t}^{1-\beta }}{4h^2}\right)\right]$$

(26)

In Equation (26), n = 1 was retained since no improvements were observed by adding more terms to the series. This means that one term is sufficient to describe the physics of this kind of test, but, in more realistic conditions, it is expected that more terms of the summation would be needed to accurately describe the solution.

5. Poromechanics Test Results and Fittings

5.1. Relaxation Tests Results—Biphasic Model

By implementing the standard procedure (i.e., confined compression, stress relaxation, no weight measurements and classic biphasic theory)—which was designed for cartilage, but it is here adopted for meniscal cylinders—the revealed mechanical parameters presented a large variability (Figure 3). Moreover, the expected trend of decreasing permeability with increasing strain was visible only for one trial (“TK11-CB-Vert” in Figure 3). This fact, coupled with high values of fitting root mean squared error (RMSE), underlines the difficulties of this approach in characterizing meniscal tissue.

5.2. Creep Test Results—Fractional Poroelastic Model

Focusing on confined compression tests, we specifically measured the decrease in weight of the cylinder during the tests and correlate this data to the amount of fluid discharged by the samples ${\mathit{j}}_f$ during the test. We now illustrate the best-fitting procedure performed between the proposed model in Equation (26) and the experimental results as shown in Figure 4a–f for a few samples taken from one of the three lateral menisci.

In particular, three free parameters are considered: $\beta ,\overline{\lambda }$, and ${\lambda }_{\beta }$. We perform the best-fitting procedure for a time window extended from t 0 to 450 s. $\gamma =0.695$ was fixed by considering specific mechanical parameters, which are summarized in

Table 4. Material parameters.

Bulk modulus	K = 1.6 × 10${}^5$ Pa
shear modulus	G = 76,923 Pa ($E=0.2\times {10}^6$ Pa, $\nu =0.3$)
Skempton coefficient	B = 0.88
Biot coefficient	$\alpha$ = 0.65
Undrained bulk modulus	$K_u$ = K/1 − $\alpha$ B

.

The values of the parameters obtained through the fittings are reported in

Table 5. The results of the fittings. The first column shows the name of the sample, the second column $\beta$ shows the order of derivative, the third column ${\lambda }_{\beta }$ is the anomalous permeability, and the fourth column specifies if the sample is extracted from a medial (M) or lateral (L) meniscus.

Sample	$\mathit{\beta }$	${\mathit{\lambda }}_{\mathit{\beta }}\left(\frac{{\mathbf{m}}^4}{\mathit{N}\times {\mathit{s}}^{1-\mathit{\beta }}}\right)$	Medial/Lateral
TK11BV	0.0242	$4.744\times {10}^{-10}$	M
TK11BR	0.0512	$5.04\times {10}^{-10}$	M
TK11BC	0.0220	$5.018\times {10}^{-10}$	M
TK16PR	0.0426	$4.619\times {10}^{-10}$	L
TK16AR	0.0259	$1.695\times {10}^{-10}$	L
TK16AV	0.0178	$0.765\times {10}^{-10}$	L
TK16PV	0.0227	$0.1083\times {10}^{-10}$	L
TK16BV	0.0397	$4.292\times {10}^{-10}$	L
TK16BC	0.0655	$10.926\times {10}^{-10}$	L
TK16BR	0.0666	$10.822\times {10}^{-10}$	L
TK17AR	0.0259	$1.331\times {10}^{-10}$	L
TK17BV	0.0553	$2.695\times {10}^{-10}$	L
TK17BC	0.0434	$4.318\times {10}^{-10}$	L
TK18BV	0.0322	$1.377\times {10}^{-10}$	M
TK18BR	0.0287	$1.897\times {10}^{-10}$	M
TK18BC	0.0289	$3.823\times {10}^{-10}$	M
TK36BC	0.0519	$7.877\times {10}^{-10}$	M
TK37BV	0.0450	$7.868\times {10}^{-10}$	L

. For samples extracted in the anterior/central/posterior and in vertical, radial and circumferential direction as shown in Figure 5.

6. Discussion

In this study, we assumed that the fluid flow is ruled by a modified version of Darcy’s law (Equation (1)). According to Equation (1), the fluid flow is not steady as modelled by the classical Darcy’s law. Instead the fluid flow rate evolves with time, more specifically it evolves with a fractional time derivative (of order $\beta$). Equation (1) also implies that the permeability is anomalous in the sense that its units are a function of the order of the derivative ${\lambda }_{\beta }=\left[\frac{L^4}{F{T}^{1-\beta }}\right]$.

The classical Darcy’s law was recovered for $\beta =0$ in Equation (1). In order to estimate how closely both models (classical and fractional Darcy’s law) fit the experimental data we compared the RMSE values. We observed in 12/18 fitted experimental tests reported in

Table 3. Complete set of the five-step test results of the central portion of the medial meniscus shown in Figure 3 for the three directions (vertical, radial and circumferential). Note that the expected trend of decreasing permeability with increasing strain was visible only for the vertical direction (TK11-CB-Vert). The root mean square error (RMSE) values are high in some cases showing that the model does not always fit the experimental behaviour.

Ramp Amplitude
Sample	Step	(% Height)	${\mathbf{H}}_{\mathbf{A}}$(MPa)	$\mathbf{k}$(${\mathbf{m}}^{\mathbf{4}}/\mathbf{s}\times \mathbf{N}$)	RMSE
	1	2	0.180	$0.858\times {10}^{-12}$	0.0003
	2	4	0.135	$2.803\times {10}^{-12}$	0.0002
TK11-CB-Circ	3	6	0.032	$0.410\times {10}^{-12}$	0.0011
	4	8	0.122	$0.953\times {10}^{-12}$	0.0004
	5	10	0.058	$1.714\times {10}^{-12}$	0.0004
	1	2	-	-	-
	2	4	0.044	$1.572\times {10}^{-12}$	0.0005
TK11-CB-Rad	3	6	0.007	$3.383\times {10}^{-12}$	0.1533
	4	8	0.037	$0.870\times {10}^{-12}$	0.0003
	5	10	0.030	$0.525\times {10}^{-12}$	0.0008
	1	2	0.0173	$1.0190\times {10}^{-12}$	0.0023
	2	4	0.0116	$0.9400\times {10}^{-12}$	0.0100
TK11-CB-Vert	3	6	0.0126	$0.8100\times {10}^{-12}$	0.0120
	4	8	0.0230	$0.4760\times {10}^{-12}$	0.0025
	5	10	0.0173	$0.4850\times {10}^{-12}$	0.0080

that the RMSE given by the fractional model was lower by 10% with respect to the classical model.

Although this might seem insufficient to justify a new theory, it is important to note that this is the first attempt to use this type of model for interpreting the rate of fluid flow inside a biological tissue. We are confident that a more precise experimental set up will give more emphasis of the benefit of having a generalized fractional model that allows for incorporation of the classical one by simply setting $\beta$ = 0.

Furthermore, it would be interesting to study the evolution of fluid flow inside the meniscus and how it varies both spatially (within the posterior/central/anterior portions) of the tissue and directionally (in the vertical/radial/circumferential directions) considering the parameters in Table 3. A first analysis of the parameters highlights that the values of the order of the fractional derivative $\beta$ and the anomalous permeability ${\lambda }_{\beta }$ are higher for the central body of the meniscus with respect of the anterior and posterior horns.

The mean values of ${\lambda }_{\beta },\beta$ in the central body are ${\lambda }_{\beta }=5.5443\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{1-\beta }}$, $\beta =0.0434$, while, in the posterior and anterior regions, are ${\lambda }_{\beta }=2.851\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{1-\beta }}$, $\beta =0.0326$ and ${\lambda }_{\beta }=1.2636\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{1-\beta }}$, $\beta =0.0232$, respectively.

Although the values of $\beta$ in the three regions might not significantly diverge from zero, they indeed affect the evolution of the fluid flow rate.

In order to investigate the role of the order of the fractional time derivative $\beta$, we specifically realized a computational simulation in which we applied a constant gradient of pressure $\parallel \nabla p\parallel$, and we calculated the fluid flow rate by applying the fractional Darcy’s law relation as given in Equation (1).

Figure 6a shows the normalized flow rate ($\overline{{\mathit{j}}_f}$) calculated considering Equation (1) in the case of a constant value of gradient of pressure $\parallel \nabla p\parallel =3\times {10}^8\frac{\mathrm{Pa}}{\mathrm{m}}$, ${\lambda }_{\beta }=5.5443\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{1-\beta }}$ and the three values of $\beta =0.0434,0.0326,0.0232$ obtained for the central body, posterior and anterior horns, respectively. Even a small value of the time fractional derivative $\beta$ affects the fluid flow response. $\beta =0$ is equivalent to considering the classical Darcy law, i.e., the fluid flow rate is constant in time.

From Figure 6a, it is possible to observe that the higher the value of $\beta$, the faster the decrease in the fluid flow rate is. Furthermore, we analysed the evolution of the fluid flow rate in the central body, posterior and anterior horns using the mean values of the parameters $\beta ,{\lambda }_{\beta }$ above. Figure 6b shows the response of the normalized flux in the three regions.

As the value of the anomalous permeability ${\lambda }_{\beta }$ in the central body is about four-times higher than in the anterior horn, the value of flux is higher in the central body and lowest in the anterior region. Moreover, as the value of $\beta$ is larger compared to the anterior and posterior horns, the fluid flux is faster (i.e., decrease in the fluid flow rate) in the central body.

We also consider how the fluid flow rate evolves in different directions, i.e., radial, vertical or circumferential. In this regard, given the paucity of data related to anterior and posterior horns, it is only possible to analyse the central body of the meniscus. The mean values of $\beta ,{\lambda }_{\beta }$ in the circumferential, radial and vertical directions are $\beta =0.0421,{\lambda }_{\beta }=6.3924\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{(1-\beta )}}$, $\beta =0.0488,{\lambda }_{\beta }=5.9196\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{(1-\beta )}}$ and $\beta =0.0393,{\lambda }_{\beta }=4.321\times {10}^{-10}\frac{{\mathrm{m}}^4}{\mathrm{N}{s}^{(1-\beta )}}$, respectively. Figure 6c pictures the evolution of the normalized fluid flow rate in the three directions. The values of $\beta$ are very close for all the three directions, differing within a few percentage points. Therefore, the decrease of flux in time (i.e., the flux velocity) can be considered as almost identical in the three main directions.

The value of $\beta$ rules the time evolution of the pore pressure diffusion in Equation (21). Increasing the value of $\beta$ implies a faster pore pressure diffusion with time. At the beginning of the test, when the fluid saturates the pores, the pore pressure carries most of the load. As the test continues, the fluid flows out of the specimen. The pore pressure decreases, and hence the solid structure starts deforming.

The rate at which the fluid flows and the pore pressure decrease and the solid structure deformation increases is ruled by the value of the fractional derivative $\beta$. A higher value of $\beta$ implies a faster pore pressure diffusion and, hence, a faster solid deformation. Figure 7 shows the evolution of the adimensional pore pressure ($\overline{p}=p(z,t)/P_A$ in Equation (21)) throughout the length of the sample along the z axis for different values of $\beta =0,0.1,0.5$ at time 15 s. It can be seen that the pore pressure drops at a faster rate with increasing $\beta$.

7. Conclusions

This work focuses on the understanding of the evolution of the fluid flow inside human meniscal tissue. We showed through $\mu$CT scans of the meniscus that the porosity significantly varies spatially within small portions of the tissue. This leads to a functionally graded permeability across the meniscal tissue. We wanted to address a question regarding how to appropriately model the fluid flow inside the meniscus. We specifically performed confined compression tests on samples extracted from three portions (posterior/central/anterior) of the structure and in three directions (vertical/radial/circumferential) in order to experimentally measure the anisotropic permeability of the human meniscal tissue.

We then correlated the weight loss of the sample with the evolution of the fluid discharged upon compression of the tissue. The results show that the weight loss of the sample is well described by a three-parameter equation derived from a fully coupled poroelastic model in which the fluid flux evolution is ruled by a generalized Darcy’s law involving fractional operators, such as derivatives of non-integer order.

We obtained the anisotropic parameters needed to characterize the fluid flow evolution inside the different parts of the meniscus, i.e., the “anomalous” permeability as well as the order of the derivative. We noticed that the parameters of this porous medium equation are functionally graded in space. Our results preliminarily suggest that the fluid flow in the central part is faster than in the posterior and anterior regions. Moreover, we noted that the flux is higher in the circumferential direction of the central body compared to the radial and vertical directions. However, the decrease in flux over time (i.e., flux velocity) can be considered almost identical in the three directions.

We believe that this work is a first attempt and a pioneering study, which experimentally investigated the anomalous behaviour of the meniscal tissue. Furthermore, the presented approach can be easily adapted to study other types of biological tissues. On the other hand, there are several limitations of this study, as here summarized:

• Additional experimental tests (unconfined compression and confined compression with control of the pressure) need to be performed in order to fully characterize the single poromechanics parameters appearing in ${\lambda }_{\beta }$ and $\overline{\lambda }$.
• A wider sample size and additional information is required to perform a reliable statistical analysis (thus far, we tested one cylinder per portion of meniscus in the vertical, radial and circumferential directions).
• This work focuses on 1D poromechanics behaviour in the three directions; however, in order to build the full anomalous permeability tensor, the coupled behaviour should be assessed in depth. Nevertheless, this study requires a more complex experimental set up.