Mathematical model of the role of degradation on matrix development in hydrogel scaffold

Abstract

Despite tremendous advances in the field of tissue engineering, a number of obstacles are still hindering its successful translation to the clinic. One of these challenges has been to design cell-laden scaffolds that can provide an appropriate environment for cells to successfully synthesize new tissue while providing a mechanical support that can resist physiological loads at the early stage of in situ implementation. A solution to this problem has been to balance tissue growth and scaffold degradation by creating new hydrogel systems that possess both hydrolytic and enzymatic degradation behaviors. Very little is known, however, about the complex behavior of these systems, emphasizing the need for a rigorous mathematical approach that can eventually assist and guide experimental advances. This paper introduces a mathematical and numerical formulation based on mixture theory, to describe the degradation, swelling, and transport of extracellular matrix (ECM) molecules released by cartilage cells (chondrocytes) within a hydrogel scaffold. The model particularly investigates the relative roles of hydrolytic and enzymatic degradations on ECM diffusion and their impacts on two important outcomes: the extent of ECM transport (and deposition) and the evolution of the scaffold’s mechanical integrity. Numerical results based on finite element show that if properly tuned, enzymatic degradation differs from hydrolytic degradation in that it can create a degradation front that is key to maintaining scaffold stiffness while allowing ECM deposition. These results therefore suggest a hydrogel design that could enable successful in situ cartilage tissue engineering.

Appendices

Appendix: Galerkin form

A mixed formulation is used with three-node elements for the solid phase and two-node elements for the other phases. The shape functions $N$ and $N_0$ and their derivatives $B,\,B_1$ and $B_0$ are defined in spherical coordinates ($R,\theta ,\psi$). The $\psi$ direction shall not be needed because of the centro-symmetric assumption (same behavior in $\theta$ and $\psi$ directions). For example, the first Piola–Kirchhoff stress $\mathbf{P}$ can be written $[P_{RR} \quad P_{\theta \theta }]^T$.

The shape functions are used to discretize the unknowns and weighting functions to derive weak form of the solid phase, alpha phase and the molecular incompressibility. They can be defined in spherical coordinates as:

$$\begin{aligned}&\mathbf{N}_l= \begin{bmatrix} N_1&N_2 \end{bmatrix} =\begin{bmatrix} \frac{1-\xi }{2}&\frac{1+\xi }{2} \end{bmatrix}\end{aligned}$$

(42)

$$\begin{aligned}&\mathbf{N}_q= \begin{bmatrix} N_1^0&N_2^0&N_3^0 \end{bmatrix} =\begin{bmatrix} \frac{\xi }{2}(\xi -1)&1-\xi ^2&\frac{\xi }{2}(\xi +1) \end{bmatrix}\end{aligned}$$

(43)

$$\begin{aligned}&\mathbf{B}_l= \begin{bmatrix} \frac{\partial N_1}{\partial R}&\frac{\partial N_2}{\partial R} \\ \frac{2 N_1}{\partial R}&\frac{2 N_2}{\partial R} \end{bmatrix}\!=\! \begin{bmatrix} \frac{-1}{l_e}&\frac{1}{l_e}\\ \frac{1-\xi }{R}&\frac{1+\xi }{R} \end{bmatrix}\end{aligned}$$

(44)

$$\begin{aligned}&\mathbf{B}_q= \begin{bmatrix} \frac{\partial N_1^0}{\partial R}&\frac{\partial N_2^0}{\partial R}&\frac{\partial N_3^0}{\partial R}\\ \frac{2 N_1^0}{ R}&\frac{2 N_2^0}{ R}&\frac{2 N_3^0}{ R} \end{bmatrix}\nonumber \\&= \begin{bmatrix} \frac{2}{l_e}(\xi -\frac{1}{2})&\frac{2}{l_e}(-2\xi )&\frac{2}{l_e}(\xi \!+\!\frac{1}{2})\\ \frac{\xi }{R}(\xi -1)&\frac{2}{R}(1-\xi ^2)&\frac{\xi }{R}(\xi \!+\!1) \end{bmatrix} \end{aligned}$$

(45)

Discretizing Eq. (27), (28) and (29) is done by injecting Eq. (30) in it. As a result, one can show that the equation of the solid (Eq. (27)) takes the following form:

$$\begin{aligned} \sum _{n=1}^{n_\mathrm{{elt}}} 4\pi \int \limits _{\Omega _0^e} R^2 \mathbf{B}_q^T \mathbf{P}^v \text{ d }R&= \sum _{n=1}^{n_\mathrm{{elt}}} 4\pi \int \limits _{\Omega _0^e} R^2 \mathbf{N}_q^T b_0 \text{ d }R \nonumber \\&\quad \!+ \sum _{n=1}^{n_\mathrm{{elt-front}}} \int \limits _{\Gamma _0^e} \mathbf{N}_q^T t_0 \text{ d }S_0 \end{aligned}$$

(46)

The alpha phase (Eq. (28)) discretized form is derived:

$$\begin{aligned}&\sum _{n=1}^{n_\mathrm{{elt}}} \left( \ \int \limits _{\Omega _0^e}\mathbf{N}_l^T \mathbf{N}_l \frac{\partial [C^\alpha ]^e}{\partial t} \text{ d }V_0 \right. + \int \limits _{\Omega _0^e}\mathbf{N}_l^T \mathbf{N}_l [\dot{\mathbf{u}}]^e \big (\mathbf{B}_l[C^\alpha ]^e_{t-\Delta t}\big )^T \text{ d }V_0 \nonumber \\&\quad \left. + \int \limits _{\Omega _0^e} k_B T \mathbf{B}_l^T \mathbf{D}^\alpha \mathbf{B}_l[C^\alpha ]^e \text{ d }V_0 + \int \limits _{\Omega _0^e} \nu ^\alpha C^\alpha \mathbf{B}_l^T \mathbf{D}^\alpha \mathbf{B}_l[\pi ]^e \text{ d }V_0 \right) \nonumber \\&\quad = \sum _{n=1}^{n_\mathrm{{elt-front}}} \left( -\int \limits _{\Gamma _0^e} \mathbf{N}_l^T I^\alpha \text{ d }S_0 \right) \end{aligned}$$

(47)

Finally, the third equation comes from the molecular incompressibility equation (Eq. (29)). Namely:

$$\begin{aligned}&\sum _{n=1}^{n_\mathrm{{elt}}} \left( \int \limits _{\Omega _0^e}(J-1) \mathbf{N}_l^T \text{ d }V_0 -\int \limits _{\Omega _0^e} \sum _\alpha \nu ^\alpha \mathbf{N}_l^T \mathbf{N}_l [C^\alpha ]^e \text{ d }V_0\right) \!=\!0\nonumber \\ \end{aligned}$$

(48)

Appendix: Linearized form

In order to solve for the unknown vector $\mathbf{y}$ in Eqs. (32), (46), (47) and (48) are linearized. Here, $\mathcal A ,\,e,\,n_\mathrm{{elt}}$ and $n_\mathrm{{elt-front}}$ denote the assembly operation, element number, number of elements and the number of elements on the boundary (frontier). The unknown vector is written as $\mathbf{y} =\begin{bmatrix} \mathbf{u}&[C^\alpha ]&[\pi ] \end{bmatrix}^T$.

1.1 Solid phase

The procedure to linearize the solid expression is to relate the first and the second Piola-Kirchhoff stresses using the expression below, and express the two terms as functions of the unknown fields (Eq (26)):

$$\begin{aligned} \mathbf{P}=\mathbf{F}\mathbf{S}\quad \text{ gives } \quad \delta \mathbf{P}=\delta \mathbf{F}\mathbf{S}+\mathbf{F}\delta \mathbf{S} \end{aligned}$$

(49)

The latter equation enables to separate the contributions from the material and geometric stiffnesses. The linearized expression thus constitutes the finite element form,

$$\begin{aligned} \mathbf{K}_s \delta \mathbf{y} = -\mathbf{R}_s \end{aligned}$$

(50)

where the stiffness $\mathbf{K}_s$ ans the residual vector $\mathbf{R}_s$ are defined as

$$\begin{aligned} \mathbf{K}_s&= \begin{bmatrix} \mathbf{K}_\mathrm{{mech}}^\mathrm{{mech}}&0&\mathbf{K}_\mathrm{{mech}}^\pi \end{bmatrix}\end{aligned}$$

(51)

$$\begin{aligned} \mathbf{R}_s&= \mathbf{F}_\mathrm{{int}} - \mathbf{F}_\mathrm{{ext}} \end{aligned}$$

(52)

The mechanical stiffness $\mathbf{K}_\mathrm{{mech}}^\mathrm{{mech}}$ is the summation of the material $\mathbf{K}_{\mathrm{{mech}},e}^\mathrm{{mat}}$ and geometrical $\mathbf{K}_{\mathrm{{mech}},e}^\mathrm{{geo}}$ contributions, respectively coming from $\delta \mathbf{S}$ and $\delta \mathbf{F}$ (Eq. (49)). $\mathbf{K}_{\mathrm{{mech}},e}^{\pi }$ comes from the osmotic pressure term in Eq. (12).

$$\begin{aligned}&\!\!\mathbf{K}_\mathrm{{mech}}^\mathrm{{mech}} = \mathcal{A }\mathbf{K}_{\mathrm{{mech}},e}^\mathrm{{mech}} = \mathcal{A } \left( \mathbf{K}_{\mathrm{{mech}},e}^\mathrm{{mat}} + \mathbf{K}_{\mathrm{{mech}},e}^\mathrm{{geo}} \right) \end{aligned}$$

(53)

$$\begin{aligned}&\!\!\mathbf{K}_\mathrm{{mech}}^\pi = \mathcal{A }\mathbf{K}_{\mathrm{{mech}},e}^\pi ;\quad \mathbf{F}_\mathrm{{int}} = \mathcal{A } \mathbf{F}_\mathrm{{int}}^e ; \quad \mathbf{F}_{\mathrm{{ext}},s} = \mathcal{A } \mathbf{F}_{\mathrm{{ext}},s}^e \end{aligned}$$

(54)

$$\begin{aligned}&\!\!\mathbf{K}_{\mathrm{{mech}},e}^\mathrm{{geo}}=4\pi \int \limits _{\Omega _0^e}\left( R^2\mathbf{B}_q^T \mathbf{S} \mathbf{B}_q \right) \text{ d }R\end{aligned}$$

(55)

$$\begin{aligned}&\!\!\mathbf{K}_{\mathrm{{mech}},e}^\mathrm{{mat}}=4\pi \int \limits _{\Omega _0^e}\left( R^2\mathbf{B}_q^T \mathbf{F} \mathbf{C} \mathbf{F} \mathbf{B}_q \right) \text{ d }R\end{aligned}$$

(56)

$$\begin{aligned}&\!\!\mathbf{K}_{\mathrm{{mech}},e}^\pi =-4\pi \int \limits _{\Omega _0^e}\left( R^2J\mathbf{B}_q^T (\mathbf{F}^{-1})^v \mathbf{N}_l \right) \text{ d }R \end{aligned}$$

(57)

$$\begin{aligned}&\!\!\mathbf{F}_\mathrm{{int}}^e = 4\pi \int \limits _{\Omega _0^e} R^2 \mathbf{B}_q^T \mathbf{P}^v \text{ d }R \end{aligned}$$

(58)

$$\begin{aligned}&\!\!\mathbf{F}_{\mathrm{{ext}},s}^e = 4\pi \int \limits _{\Omega _0^e} R^2 \mathbf{N}_q^T b_0 \text{ d }R + \int \limits _{\Gamma _0^e} \mathbf{N}_q^T t_0 \text{ d }S_0 \end{aligned}$$

(59)

1.2 Alpha phase

The alpha phase equation is linearized after discretization. Thus, Eq. (47) gives

$$\begin{aligned} \mathbf{C}_\alpha \dot{\mathbf{y}} + \mathbf{K}_\alpha \delta \mathbf{y} = -\mathbf{R}_\alpha \end{aligned}$$

(60)

where the damping matrix $\mathbf{C}_\alpha$, stiffness matrix $\mathbf{K}_\alpha$ and residual vector $\mathbf{R}_\alpha$ are defined as follows

$$\begin{aligned} \mathbf{C}_\alpha&= \begin{bmatrix} \mathbf{C}_{\alpha }^\mathrm{{mech}}&\mathbf{C}_\alpha ^\alpha&0 \end{bmatrix}\end{aligned}$$

(61)

$$\begin{aligned} \mathbf{K}_\alpha&= \begin{bmatrix} 0&\mathbf{K}_\alpha ^\alpha&\mathbf{K}_\alpha ^\pi \end{bmatrix}\end{aligned}$$

(62)

$$\begin{aligned} \mathbf{R}_\alpha&= \mathbf{F}_\mathrm{{int},\alpha } - \mathbf{F}_\mathrm{{ext},\alpha } \end{aligned}$$

(63)

The damping stiffnesses ($\mathbf{C}_\alpha ^\mathrm{{mech}}$ and $\mathbf{C}_\alpha ^\alpha$) come from the time dependency in Eq. (15), and the stiffness matrices $\mathbf{K}_\alpha ^\alpha$ and $\mathbf{K}_\alpha ^\pi$ represent the stiffnesses from Ficks’s and Darcy’s law (Eq. 20).

$$\begin{aligned}&\!\!\mathbf{C}_{\alpha }^{\alpha } = \mathcal{A } \mathbf{C}_{\alpha ,e}^{\alpha }; \quad \mathbf{C}_{\alpha }^\mathrm{{mech}} = \mathcal{A } \mathbf{C}_{\alpha ,e}^\mathrm{{mech}} ;\quad \mathbf{K}_{\alpha }^{\alpha } = \mathcal A \mathbf{K}_{\alpha ,e}^{\alpha }\end{aligned}$$

(64)

$$\begin{aligned}&\!\!\mathbf{F}_\mathrm{{ext},\alpha } = \mathcal{A } \mathbf{F}_\mathrm{{ext},\alpha }^e ;\quad \mathbf{F}_\mathrm{{int},\alpha } = \mathcal A \mathbf{F}_\mathrm{{int},\alpha }^e \end{aligned}$$

(65)

$$\begin{aligned}&\!\!\mathbf{K}_{\alpha ,e}^{\alpha } = 4\pi \int \limits _{\Omega _0^e} k_B T \mathbf{B}_l^T \mathbf{D}^\alpha \mathbf{B}_l \text{ d }R \end{aligned}$$

(66)

$$\begin{aligned}&\!\!\mathbf{K}_{\alpha ,e}^{\pi }= 4\pi \int \limits _{\Omega _0^e} \nu ^\alpha C^\alpha \mathbf{B}_l^T \mathbf{D}^\alpha \mathbf{B}_l \text{ d }R\end{aligned}$$

(67)

$$\begin{aligned}&\!\!\mathbf{C}_{\alpha ,e}^{\alpha } = 4\pi \int \limits _{\Omega _0^e}\mathbf{N}_l^T \mathbf{N}_l \text{ d }R\end{aligned}$$

(68)

$$\begin{aligned}&\!\!\mathbf{C}_{\alpha ,e}^\mathrm{{mech}} = 4\pi \int \limits _{\Omega _0^e}\mathbf{N}_l^T \mathbf{N}_l (\mathbf{B}_l[C^\alpha ]^e_{t-\Delta t})^T \text{ d }R\end{aligned}$$

(69)

$$\begin{aligned}&\!\!\mathbf{F}_\mathrm{{ext},\alpha }^e = -\int \limits _{\Gamma _0^e} \mathbf{N}_l^T I^\alpha \text{ d }S_0\end{aligned}$$

(70)

$$\begin{aligned}&\!\!\mathbf{F}_\mathrm{{int},\alpha }^e = \left( 4\pi \int \limits _{\Omega _0^e} k_B T \mathbf{B}_l^T \mathbf{D}^\alpha \mathbf{B}_l \text{ d }R \right) C^\alpha \nonumber \\&\qquad \qquad + \left( 4\pi \int \limits _{\Omega _0^e} \nu ^\alpha C^\alpha \mathbf{B}_l^T \mathbf{D}^\alpha \mathbf{B}_l \right) \pi \end{aligned}$$

(71)

1.3 Molecular incompressibility

Finally, Eq. (48) linearization results in the last equation:

$$\begin{aligned} \mathbf{K}_\mathrm{{constr}} \delta \mathbf{y} = -\mathbf{R}_\mathrm{{constr}} \end{aligned}$$

(72)

where the stiffness matrix $\mathbf{K}_\mathrm{{constr}}$ and residual vector $\mathbf{R}_\mathrm{{constr}}$ are defined as follows

$$\begin{aligned} \mathbf{K}_\mathrm{{constr}}&= \begin{bmatrix} \mathbf{K}_\mathrm{{constr}}^\mathrm{{mech}}&\mathbf{K}_\mathrm{{constr}}^\alpha&0 \end{bmatrix}\end{aligned}$$

(73)

$$\begin{aligned} \mathbf{R}_\mathrm{{constr}}&= \mathbf{F}_\mathrm{{int,constr}} - \mathbf{F}_\mathrm{{ext,constr}} \end{aligned}$$

(74)

The stiffness matrix $\mathbf{K}_\mathrm{{constr}}$ emphasizes that the constraint has an impact on both the solid ($\mathbf{K}_\mathrm{{constr}}^\mathrm{{mech}}$) and alpha phase ($\mathbf{K}_\mathrm{{constr}}^\alpha$). It comes, respectively, from the term $J$ and $\nu ^\alpha C^\alpha$ in Eq. (9).

$$\begin{aligned}&\mathbf{K}_\mathrm{{constr}}^\mathrm{{mech}} = \mathcal{A }\mathbf{K}_{\mathrm{{constr}},e}^\mathrm{{mech}}\end{aligned}$$

(75)

$$\begin{aligned}&\mathbf{K}_\mathrm{{constr}}^\alpha = \mathcal{A }\mathbf{K}_{\mathrm{{constr}},e}^\alpha \end{aligned}$$

(76)

$$\begin{aligned}&\mathbf{K}_{\mathrm{{constr}},e}^\mathrm{{mech}}=4\pi \int \limits _{\Omega _0^e} \left( J \mathbf{N}_l^T [(\mathbf{F}^{-T})^v]^T \mathbf{B}_q \right) \text{ d }R\end{aligned}$$

(77)

$$\begin{aligned}&\mathbf{K}_{\mathrm{{constr}},e}^\alpha =-4\pi \int \limits _{\Omega _0^e}\left( \sum _\alpha \nu ^\alpha \mathbf{N}_l^T \mathbf{N}_l\right) \text{ d }R\end{aligned}$$

(78)

$$\begin{aligned}&\mathbf{F}^e_\mathrm{{int,constr}} = 4\pi \int \limits _{\Omega _0^e}\left( (J-1) \mathbf{N}_l^T \right) \text{ d }R \nonumber \\&\qquad \qquad \qquad - 4\pi \int \limits _{\Omega _0^e} \left( \sum _\alpha \nu ^\alpha \mathbf{N}_l^T \mathbf{N}_l [C^\alpha ]^e \right) \text{ d }R\end{aligned}$$

(79)

$$\begin{aligned}&\mathbf{F}^e_\mathrm{{ext,constr}} = 0 \end{aligned}$$

(80)