Mathematical Modeling of Biological Tissue Regeneration
The epidermal wound healing process can be modeled using a system of reaction-diffusion equations that describe the spatiotemporal dynamics of cell density, growth factors, and extracellular matrix components.
When \( p = 0 \) and \( s_c = 0 \), the PDE reduces to the heat equation:
$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial r^2}$$We decompose \( u(r,t) = u_s(r) + v(r,t) \), where \( u_s(r) \) is the steady-state solution satisfying:
$$D \frac{d^2 u_s}{dr^2} = 0, \quad u_s(0) = 1, \quad \frac{du_s}{dr}(r_0) = 0$$Solving this gives \( u_s(r) = 1 \). Thus \( v(r,t) = u(r,t) - 1 \) satisfies:
$$\frac{\partial v}{\partial t} = D \frac{\partial^2 v}{\partial r^2}, \quad v(0,t) = 0, \quad \frac{\partial v}{\partial r}(r_0,t) = 0, \quad v(r,0) = -1$$Using separation of variables \( v(r,t) = R(r)T(t) \):
$$\frac{T'}{DT} = \frac{R''}{R} = -\lambda^2$$The spatial ODE with boundary conditions yields:
$$R_m(r) = \sin\left( \frac{(2m-1)\pi r}{2r_0} \right), \quad \lambda_m = \frac{(2m-1)\pi}{2r_0}, \quad m = 1, 2, \dots$$ $$T_m(t) = e^{-D \lambda_m^2 t}$$Thus the general solution becomes:
$$v(r,t) = \sum_{m=1}^{\infty} c_m e^{-D\lambda_m^2 t} \sin\left( \frac{(2m-1)\pi r}{2r_0} \right)$$Applying the initial condition \( v(r,0) = -1 \):
$$-1 = \sum_{m=1}^{\infty} c_m \sin\left( \frac{(2m-1)\pi r}{2r_0} \right)$$ $$c_m = -\frac{4}{(2m-1)\pi}$$Finally, the complete analytical solution is:
$$\boxed{u(r,t) = 1 - \sum_{m=1}^{\infty} \frac{4}{(2m-1)\pi} \exp\left(-D \left( \frac{(2m-1)\pi}{2r_0} \right)^2 t \right) \sin\left( \frac{(2m-1)\pi r}{2r_0} \right)}$$The expanded PDE is discretized as:
\[ \frac{\partial u}{\partial t} \approx \frac{u_i^{n+1} - u_i^n}{\Delta t} \] \[ \frac{\partial^2 u}{\partial r^2} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta r)^2}, \quad \left( \frac{\partial u}{\partial r} \right)^2 \approx \left( \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta r} \right)^2 \]
The update equation for interior points \(i = 2, \dots, N-1\) is:
\[ \begin{aligned} u_i^{n+1} = u_i^n + \Delta t \bigg[ D \bigg( \underbrace{\left(1 - \frac{u_i^n}{u_0}\right)^p \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta r)^2}}_{\text{Diffusion}} - \underbrace{\frac{p}{u_0} \left(1 - \frac{u_i^n}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^n - u_{i-1}^n}{2 \Delta r} \right)^2}_{\text{Nonlinear flux}} \bigg) + \underbrace{s_c u_i^n \left(1 - \frac{u_i^n}{u_0}\right)}_{\text{Source}} \bigg] \end{aligned} \]
Boundary conditions:
\[ \begin{aligned} u_1^{n+1} &= 1 \quad \text{(Dirichlet)} \\ u_N^{n+1} &= u_{N-1}^{n+1} \quad \text{(Neumann)} \end{aligned} \]
Stability requires:
\[ \Delta t \leq \frac{(\Delta r)^2}{2D} \]
Discretize at \( t_{n+1} \):
\[ \begin{aligned} \frac{u_i^{n+1} - u_i^n}{\Delta t} = D \bigg[ &\left(1 - \frac{u_i^{n+1}}{u_0}\right)^p \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{(\Delta r)^2} \\ &- \frac{p}{u_0} \left(1 - \frac{u_i^{n+1}}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^{n+1} - u_{i-1}^{n+1}}{2 \Delta r} \right)^2 \bigg] \\ &+ s_c u_i^{n+1} \left(1 - \frac{u_i^{n+1}}{u_0}\right) \end{aligned} \]
This forms a nonlinear system:
\[ F_i(\mathbf{u}^{n+1}) = u_i^{n+1} - u_i^n - \Delta t \cdot \text{RHS}(\mathbf{u}^{n+1}) = 0 \]
Solved with Newton-Raphson:
\[ J(\mathbf{u}_k) \delta \mathbf{u} = -\mathbf{F}(\mathbf{u}_k), \quad \mathbf{u}_{k+1} = \mathbf{u}_k + \delta \mathbf{u} \]
where \( J_{ij} = \partial F_i/\partial u_j^{n+1} \) is the Jacobian.
\[ \begin{aligned} \frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{1}{2} \Big[ & D \left( \left(1 - \frac{u_i^n}{u_0}\right)^p \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta r)^2} - \frac{p}{u_0} \left(1 - \frac{u_i^n}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta r} \right)^2 \right) \\ + & s_c u_i^n \left(1 - \frac{u_i^n}{u_0}\right) \Big] \\ + \frac{1}{2} \Big[ & D \left( \left(1 - \frac{u_i^{n+1}}{u_0}\right)^p \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{(\Delta r)^2} - \frac{p}{u_0} \left(1 - \frac{u_i^{n+1}}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^{n+1} - u_{i-1}^{n+1}}{2\Delta r} \right)^2 \right) \\ + & s_c u_i^{n+1} \left(1 - \frac{u_i^{n+1}}{u_0}\right) \Big] \end{aligned} \]
Solved as a nonlinear system using Newton-Raphson. Second-order accurate and unconditionally stable.The Method of Lines discretizes spatial derivatives to convert the PDE into a system of ODEs.
Define grid points \( r_i = (i-1)\Delta r \) with \( \Delta r = \frac{r_0}{N-1} \). Let \( u_i(t) = u(r_i, t) \).
For interior points (\( i=2,\dots,N-1 \)):
\[ \frac{du_i}{dt} = D \left[ \left(1 - \frac{u_i}{u_0}\right)^p \frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta r)^2} - \frac{p}{u_0} \left(1 - \frac{u_i}{u_0}\right)^{p-1} \left( \frac{u_{i+1} - u_{i-1}}{2\Delta r} \right)^2 \right] + s_c u_i \left(1 - \frac{u_i}{u_0}\right) \]
\[ u_1(t) = 1 \implies \frac{du_1}{dt} = 0 \]
\[ \frac{u_{N+1} - u_{N-1}}{2\Delta r} = 0 \implies u_{N+1} = u_{N-1} \]
Then: \[ \frac{du_N}{dt} = D \left(1 - \frac{u_N}{u_0}\right)^p \frac{2(u_{N-1} - u_N)}{(\Delta r)^2} + s_c u_N \left(1 - \frac{u_N}{u_0}\right) \]
\[ u_i(0) = \begin{cases} 1 & i=1 \\ 0 & i=2,\dots,N \end{cases} \]
Solve the ODE system:
\[ \frac{d\mathbf{u}}{dt} = \mathbf{F}(t, \mathbf{u}), \quad \mathbf{u}(0) = \mathbf{u}_0 \]
using a stiff ODE solver (e.g., BDF method).
Multiply by test function \( v(r) \) with \( v(0) = 0 \):
\[ \int_0^{r_0} \frac{\partial u}{\partial t} v \, dr = \int_0^{r_0} \left[ D \frac{\partial}{\partial r} \left( \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \right) + s_c u \left(1 - \frac{u}{u_0}\right) \right] v \, dr \]
Integrate diffusion term by parts:
\[ \int_0^{r_0} D \frac{\partial}{\partial r} \left( \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \right) v \, dr = \left[ D \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} v \right]_0^{r_0} - \int_0^{r_0} D \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \frac{\partial v}{\partial r} \, dr \]
Boundary term vanishes due to \( v(0) = 0 \) and \( \frac{\partial u}{\partial r}(r_0) = 0 \), yielding:
\[ \boxed{ \int_0^{r_0} \frac{\partial u}{\partial t} v \, dr + \int_0^{r_0} D \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \frac{\partial v}{\partial r} \, dr = \int_0^{r_0} s_c u \left(1 - \frac{u}{u_0}\right) v \, dr } \]
Discretize domain \( [0, r_0] \) into \( N \) nodes with piecewise linear basis functions \( \phi_j(r) \). Approximate solution:
\[ u_h(r, t) = \phi_1(r) + \sum_{j=2}^N u_j(t) \phi_j(r) \]
where \( u_1(t) = 1 \) (Dirichlet condition).
For test functions \( v = \phi_i \) (\( i = 2, \dots, N \)):
\[ \sum_{j=2}^N \frac{du_j}{dt} \underbrace{\int_0^{r_0} \phi_i \phi_j \, dr}_{M_{ij}} + \sum_{j=2}^N u_j(t) \underbrace{\int_0^{r_0} D \left(1 - \frac{u_h}{u_0}\right)^p \frac{d\phi_i}{dr} \frac{d\phi_j}{dr} \, dr}_{K_{ij}(u)} = \underbrace{\int_0^{r_0} s_c u_h \left(1 - \frac{u_h}{u_0}\right) \phi_i \, dr}_{f_i(u)} \]
Resulting in:
\[ M \frac{d\mathbf{u}}{dt} + K(\mathbf{u}) \mathbf{u} = \mathbf{f}(\mathbf{u}) \]
where \( \mathbf{u} = [u_2, \dots, u_N]^T \).
\[ M \frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\Delta t} + K(\mathbf{u}^{n+1}) \mathbf{u}^{n+1} = \mathbf{f}(\mathbf{u}^{n+1}) \]
Rearrange to nonlinear system:
\[ \mathbf{G}(\mathbf{u}^{n+1}) = M \mathbf{u}^{n+1} + \Delta t K(\mathbf{u}^{n+1}) \mathbf{u}^{n+1} - \Delta t \mathbf{f}(\mathbf{u}^{n+1}) - M \mathbf{u}^n = \mathbf{0} \]
Solved with Newton-Raphson method.
Explore a detailed simulation of the model, showcasing its behavior under various conditions.
View SimulationWe are a team of passionate students and engineers dedicated to applying mathematical modeling, simulation, and numerical techniques to solve real-world biomedical and engineering challenges.
Team Members:
The epidermal wound healing process can be modeled using a system of reaction-diffusion equations that describe the spatiotemporal dynamics of cell density, growth factors, and extracellular matrix components.
When \( p = 0 \) and \( s_c = 0 \), the PDE reduces to the heat equation:
$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial r^2}$$We decompose \( u(r,t) = u_s(r) + v(r,t) \), where \( u_s(r) \) is the steady-state solution satisfying:
$$D \frac{d^2 u_s}{dr^2} = 0, \quad u_s(0) = 1, \quad \frac{du_s}{dr}(r_0) = 0$$Solving this gives \( u_s(r) = 1 \). Thus \( v(r,t) = u(r,t) - 1 \) satisfies:
$$\frac{\partial v}{\partial t} = D \frac{\partial^2 v}{\partial r^2}, \quad v(0,t) = 0, \quad \frac{\partial v}{\partial r}(r_0,t) = 0, \quad v(r,0) = -1$$Using separation of variables \( v(r,t) = R(r)T(t) \):
$$\frac{T'}{DT} = \frac{R''}{R} = -\lambda^2$$The spatial ODE with boundary conditions yields:
$$R_m(r) = \sin\left( \frac{(2m-1)\pi r}{2r_0} \right), \quad \lambda_m = \frac{(2m-1)\pi}{2r_0}, \quad m = 1, 2, \dots$$ $$T_m(t) = e^{-D \lambda_m^2 t}$$Thus the general solution becomes:
$$v(r,t) = \sum_{m=1}^{\infty} c_m e^{-D\lambda_m^2 t} \sin\left( \frac{(2m-1)\pi r}{2r_0} \right)$$Applying the initial condition \( v(r,0) = -1 \):
$$-1 = \sum_{m=1}^{\infty} c_m \sin\left( \frac{(2m-1)\pi r}{2r_0} \right)$$ $$c_m = -\frac{4}{(2m-1)\pi}$$Finally, the complete analytical solution is:
$$\boxed{u(r,t) = 1 - \sum_{m=1}^{\infty} \frac{4}{(2m-1)\pi} \exp\left(-D \left( \frac{(2m-1)\pi}{2r_0} \right)^2 t \right) \sin\left( \frac{(2m-1)\pi r}{2r_0} \right)}$$The expanded PDE is discretized as:
\[ \frac{\partial u}{\partial t} \approx \frac{u_i^{n+1} - u_i^n}{\Delta t} \] \[ \frac{\partial^2 u}{\partial r^2} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta r)^2}, \quad \left( \frac{\partial u}{\partial r} \right)^2 \approx \left( \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta r} \right)^2 \]
The update equation for interior points \(i = 2, \dots, N-1\) is:
\[ \begin{aligned} u_i^{n+1} = u_i^n + \Delta t \bigg[ D \bigg( \underbrace{\left(1 - \frac{u_i^n}{u_0}\right)^p \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta r)^2}}_{\text{Diffusion}} - \underbrace{\frac{p}{u_0} \left(1 - \frac{u_i^n}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^n - u_{i-1}^n}{2 \Delta r} \right)^2}_{\text{Nonlinear flux}} \bigg) + \underbrace{s_c u_i^n \left(1 - \frac{u_i^n}{u_0}\right)}_{\text{Source}} \bigg] \end{aligned} \]
Boundary conditions:
\[ \begin{aligned} u_1^{n+1} &= 1 \quad \text{(Dirichlet)} \\ u_N^{n+1} &= u_{N-1}^{n+1} \quad \text{(Neumann)} \end{aligned} \]
Stability requires:
\[ \Delta t \leq \frac{(\Delta r)^2}{2D} \]
Discretize at \( t_{n+1} \):
\[ \begin{aligned} \frac{u_i^{n+1} - u_i^n}{\Delta t} = D \bigg[ &\left(1 - \frac{u_i^{n+1}}{u_0}\right)^p \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{(\Delta r)^2} \\ &- \frac{p}{u_0} \left(1 - \frac{u_i^{n+1}}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^{n+1} - u_{i-1}^{n+1}}{2 \Delta r} \right)^2 \bigg] \\ &+ s_c u_i^{n+1} \left(1 - \frac{u_i^{n+1}}{u_0}\right) \end{aligned} \]
This forms a nonlinear system:
\[ F_i(\mathbf{u}^{n+1}) = u_i^{n+1} - u_i^n - \Delta t \cdot \text{RHS}(\mathbf{u}^{n+1}) = 0 \]
Solved with Newton-Raphson:
\[ J(\mathbf{u}_k) \delta \mathbf{u} = -\mathbf{F}(\mathbf{u}_k), \quad \mathbf{u}_{k+1} = \mathbf{u}_k + \delta \mathbf{u} \]
where \( J_{ij} = \partial F_i/\partial u_j^{n+1} \) is the Jacobian.
\[ \begin{aligned} \frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{1}{2} \Big[ & D \left( \left(1 - \frac{u_i^n}{u_0}\right)^p \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta r)^2} - \frac{p}{u_0} \left(1 - \frac{u_i^n}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta r} \right)^2 \right) \\ + & s_c u_i^n \left(1 - \frac{u_i^n}{u_0}\right) \Big] \\ + \frac{1}{2} \Big[ & D \left( \left(1 - \frac{u_i^{n+1}}{u_0}\right)^p \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{(\Delta r)^2} - \frac{p}{u_0} \left(1 - \frac{u_i^{n+1}}{u_0}\right)^{p-1} \left( \frac{u_{i+1}^{n+1} - u_{i-1}^{n+1}}{2\Delta r} \right)^2 \right) \\ + & s_c u_i^{n+1} \left(1 - \frac{u_i^{n+1}}{u_0}\right) \Big] \end{aligned} \]
Solved as a nonlinear system using Newton-Raphson. Second-order accurate and unconditionally stable.The Method of Lines discretizes spatial derivatives to convert the PDE into a system of ODEs.
Define grid points \( r_i = (i-1)\Delta r \) with \( \Delta r = \frac{r_0}{N-1} \). Let \( u_i(t) = u(r_i, t) \).
For interior points (\( i=2,\dots,N-1 \)):
\[ \frac{du_i}{dt} = D \left[ \left(1 - \frac{u_i}{u_0}\right)^p \frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta r)^2} - \frac{p}{u_0} \left(1 - \frac{u_i}{u_0}\right)^{p-1} \left( \frac{u_{i+1} - u_{i-1}}{2\Delta r} \right)^2 \right] + s_c u_i \left(1 - \frac{u_i}{u_0}\right) \]
\[ u_1(t) = 1 \implies \frac{du_1}{dt} = 0 \]
\[ \frac{u_{N+1} - u_{N-1}}{2\Delta r} = 0 \implies u_{N+1} = u_{N-1} \]
Then: \[ \frac{du_N}{dt} = D \left(1 - \frac{u_N}{u_0}\right)^p \frac{2(u_{N-1} - u_N)}{(\Delta r)^2} + s_c u_N \left(1 - \frac{u_N}{u_0}\right) \]
\[ u_i(0) = \begin{cases} 1 & i=1 \\ 0 & i=2,\dots,N \end{cases} \]
Solve the ODE system:
\[ \frac{d\mathbf{u}}{dt} = \mathbf{F}(t, \mathbf{u}), \quad \mathbf{u}(0) = \mathbf{u}_0 \]
using a stiff ODE solver (e.g., BDF method).
Multiply by test function \( v(r) \) with \( v(0) = 0 \):
\[ \int_0^{r_0} \frac{\partial u}{\partial t} v \, dr = \int_0^{r_0} \left[ D \frac{\partial}{\partial r} \left( \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \right) + s_c u \left(1 - \frac{u}{u_0}\right) \right] v \, dr \]
Integrate diffusion term by parts:
\[ \int_0^{r_0} D \frac{\partial}{\partial r} \left( \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \right) v \, dr = \left[ D \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} v \right]_0^{r_0} - \int_0^{r_0} D \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \frac{\partial v}{\partial r} \, dr \]
Boundary term vanishes due to \( v(0) = 0 \) and \( \frac{\partial u}{\partial r}(r_0) = 0 \), yielding:
\[ \boxed{ \int_0^{r_0} \frac{\partial u}{\partial t} v \, dr + \int_0^{r_0} D \left(1 - \frac{u}{u_0}\right)^p \frac{\partial u}{\partial r} \frac{\partial v}{\partial r} \, dr = \int_0^{r_0} s_c u \left(1 - \frac{u}{u_0}\right) v \, dr } \]
Discretize domain \( [0, r_0] \) into \( N \) nodes with piecewise linear basis functions \( \phi_j(r) \). Approximate solution:
\[ u_h(r, t) = \phi_1(r) + \sum_{j=2}^N u_j(t) \phi_j(r) \]
where \( u_1(t) = 1 \) (Dirichlet condition).
For test functions \( v = \phi_i \) (\( i = 2, \dots, N \)):
\[ \sum_{j=2}^N \frac{du_j}{dt} \underbrace{\int_0^{r_0} \phi_i \phi_j \, dr}_{M_{ij}} + \sum_{j=2}^N u_j(t) \underbrace{\int_0^{r_0} D \left(1 - \frac{u_h}{u_0}\right)^p \frac{d\phi_i}{dr} \frac{d\phi_j}{dr} \, dr}_{K_{ij}(u)} = \underbrace{\int_0^{r_0} s_c u_h \left(1 - \frac{u_h}{u_0}\right) \phi_i \, dr}_{f_i(u)} \]
Resulting in:
\[ M \frac{d\mathbf{u}}{dt} + K(\mathbf{u}) \mathbf{u} = \mathbf{f}(\mathbf{u}) \]
where \( \mathbf{u} = [u_2, \dots, u_N]^T \).
\[ M \frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\Delta t} + K(\mathbf{u}^{n+1}) \mathbf{u}^{n+1} = \mathbf{f}(\mathbf{u}^{n+1}) \]
Rearrange to nonlinear system:
\[ \mathbf{G}(\mathbf{u}^{n+1}) = M \mathbf{u}^{n+1} + \Delta t K(\mathbf{u}^{n+1}) \mathbf{u}^{n+1} - \Delta t \mathbf{f}(\mathbf{u}^{n+1}) - M \mathbf{u}^n = \mathbf{0} \]
Solved with Newton-Raphson method.
Explore a detailed simulation of the model, showcasing its behavior under various conditions.
We are a team of passionate students and engineers dedicated to applying mathematical modeling, simulation, and numerical techniques to solve real-world biomedical and engineering challenges.