Reaction-Diffusion Model

Introduction

Diabetic chronic wounds are one of the most severe and challenging complications of diabetes, affecting up to 25% of patients and characterized by a persistent failure to heal that often leads to infection, amputation, and even death^[1]. The pathophysiology underlying the difficulty in healing these wounds is extremely complex, with two core obstacles being a persistent chronic inflammatory state and severe impairment of angiogenesis^[2]. The former prevents the wound from transitioning from the inflammatory to the reparative phase, while the latter leads to a lack of necessary oxygen and nutrient supply to the nascent tissue.

To address both challenges simultaneously, we propose an innovative "living therapy" strategy: encapsulating genetically engineered yeast cells capable of synergistically secreting the anti-inflammatory factor IL-4 and the pro-angiogenic factor VEGF, along with antimicrobial peptides, in a biocompatible chitosan hydrogel. This system acts as an in-situ "micro-bioreactor," designed to achieve a sustained and stable release of therapeutic proteins directly at the wound site. In the hostile diabetic microenvironment, the body's own healing signals are suppressed and topically applied growth factors are rapidly degraded due to high levels of proteases like matrix metalloproteinases (MMPs)^[3]. By providing a continuous, high-concentration local supply, our system is designed to overwhelm these inhibitory effects and re-establish the conditions necessary for repair, a significant advantage over single-dose administrations that struggle with short half-lives and the need for frequent dressing changes.

However, moving from theoretical design to clinical application, a core scientific question must be answered: In the harsh microenvironment specific to diabetic chronic wounds, how long does it take for the yeast-secreted IL-4 and VEGF to reach and maintain effective therapeutic concentrations within the granulation tissue that genuinely needs repair? The answer to this question is crucial as it directly determines the feasibility and clinical potential of this therapy. Diabetic chronic wounds exudate contains high concentrations of matrix metalloproteinases (MMPs), which can efficiently degrade various growth factors, including VEGF^[3].

Therefore, we employ computational modeling to construct a reaction-diffusion partial differential equation (PDE) model that can accurately describe the dynamic changes of substances in space and time. PDEs are the gold standard language for describing physicochemical processes such as diffusion, reaction, and convection. By creating a "digital twin" of the wound microenvironment, we are able to:

1. 1. Simulate complex biochemical processes in a computer to quantify the impact of adverse factors like MMPs.
2. 2. Conduct rapid "virtual experiments" to predict the "onset time" of the drug.
3. 3. Identify key bottleneck parameters affecting efficacy, thereby providing clear and focused validation targets for subsequent wet-lab experiments.

The purpose of this study is precisely to combine our wet-lab experimental goal (i.e., verifying the efficacy of the yeast-hydrogel system) with rigorous theoretical predictions by establishing and solving such a high-fidelity PDE model, enabling efficient experimental research and development guided by computational simulation.

Model Assumptions

To translate the complex biological problem into a solvable mathematical model, we make the following key assumptions and justify their reasonableness:

1. 1. One-Dimensional Geometry Simplification (1D Geometry):
   • Assumption: We simplify the wound into a one-dimensional spatial domain, which represents the depth profile from the outer surface of the hydrogel (x = -H) vertically down to the tissue (x = L). Since our goal is to treat deeper diabetic wounds (reaching the granulation layer), we can set the interface at x = 0 to be in direct contact with the granulation tissue.
   • Rationale: The main concentration gradient of drug delivery exists in the vertical direction. Combined with the clinical scenario of interest in this study (deep or debrided diabetic wounds where the dressing directly covers the granulation layer), drug exposure in the near-interface region is most clinically relevant.
2. 2. Continuum Hypothesis:
   • Assumption: We treat molecules like IL-4 and VEGF as continuous concentration fields \(C(x,t)\), rather than discrete individual molecules.
   • Rationale: On the micrometer to millimeter scale we are concerned with, the number of molecules is extremely large, and their collective behavior can be accurately described by a macroscopic average concentration. This is the foundation of all reaction-diffusion models.
3. 3. Constant Yeast Production:
   • Assumption: We assume that the yeast encapsulated in the hydrogel will produce therapeutic factors at a constant average rate \(P_i\) after a brief adaptation period.
   • Rationale: This is an initial simplification of the complex yeast physiology. It allows us to focus on the drug delivery process itself, rather than the growth dynamics of the yeast. This production rate \(P_i\) is a key design parameter that can be precisely measured in subsequent experiments and used as an input for the model.
4. 4. Quasi-Steady-State Assumption:
   • Assumption: The binding and dissociation reactions of therapeutic factors with the hydrogel matrix and albumin are much faster than their diffusion in space or their degradation.
   • Rationale: Intermolecular binding/dissociation typically occurs on the nanosecond to millisecond scale, whereas diffusion across hundreds of micrometers takes hours. QSSA allows us to simplify these rapid intermediate steps into algebraic equilibrium relationships, thus avoiding the need to solve additional, potentially numerically unstable ordinary differential equations. This makes the model easier to solve while retaining its core physical meaning.
5. 5. Homogenized Tissue Layers:
   • Assumption: In this scenario, the interface is the surface of the granulation layer (the thickness of the necrotic layer is approximately 0), and the tissue domain is divided into two layers: the granulation layer and healthy tissue. The diffusion coefficient \(D\) and the effective consumption rate \(k_{eff}\) are piecewise constant within each layer.
   • Rationale: This setting reflects the clinical reality of "a deep wound having reached the granulation layer," focusing the criterion on the time to reach the target concentration in the near-interface region.
6. 6. Constant Exudate Source:
   • Assumption: The concentrations of MMPs and albumin at the wound interface remain constant over the simulated time frame.
   • Rationale: Chronic wounds continuously produce a large amount of exudate, forming a huge "buffer pool." We therefore assume that the small amount of MMPs and albumin diffusing from this source into the hydrogel will not significantly alter the source concentration.

Together, these assumptions construct an idealized model that both reflects key biological realities and is mathematically tractable, making it a powerful tool for hypothesis testing and parameter sensitivity analysis.

Model Development

To quantitatively predict the spatiotemporal dynamics of therapeutic factors from the engineered yeast-hydrogel system within the diabetic wound microenvironment, we established a system of coupled reaction-diffusion partial differential equations (PDEs) based on the law of mass conservation. The model simplifies the complex wound environment into a one-dimensional (1D) computational domain, representing the depth profile extending vertically from the hydrogel's outer surface at \(x = -H\) down to the tissue depth at \(x = L\), with the gel-tissue interface located at \(x = 0\). This geometry effectively captures the primary mass transport gradients perpendicular to the wound surface. We simultaneously solve for the concentration distributions of four key species: IL-4 \(C_{IL4}(x,t)\), VEGF \(C_{VEGF}(x,t)\), matrix metalloproteinases \(C_{M}(x,t)\), and albumin \(C_{A}(x,t)\).

Dynamic Model within the Hydrogel Domain (Ωg)

The hydrogel domain, $-H \le x < 0$, serves as the "living bioreactor," where the mass transport and reaction network are most complex. First, MMPs and albumin from the wound exudate infiltrate from the boundary at \(x=0\). Their diffusion and reaction processes are described by the following equations:

$$ \frac{\partial C_M}{\partial t} = D_M \frac{\partial^2 C_M}{\partial x^2} - k_{M,inact} C_M \quad $$ $$ \frac{\partial C_A}{\partial t} = D_A \frac{\partial^2 C_A}{\partial x^2} \quad $$

Here, MMPs undergo first-order inactivation (with rate constant \(k_{M,inact}\)), while albumin is considered an inert diffusing species in this domain.

For the therapeutic factors IL-4 and VEGF, their fate within the hydrogel is governed by production, diffusion, degradation, and reversible binding to both the hydrogel matrix and the infiltrating albumin. Taking IL-4 as an example, the complete system of dynamic equations for its free concentration \(C_{IL4}\), chitosan-bound concentration \(B_{IL4}^{chito}\), and albumin-bound concentration \(B_{IL4}^{alb}\) can be written as:

$$ \begin{aligned} \frac{\partial C_{IL4}}{\partial t} = & D_{IL4,g} \frac{\partial^2 C_{IL4}}{\partial x^2} + P_{IL4} - (k_{IL4,g}^{deg} + k_{IL4,M} C_M) C_{IL4} \\ & - \left( k_{chito,IL4}^+ C_{IL4} (S_{chito} - B_{IL4}^{chito}) - k_{chito,IL4}^- B_{IL4}^{chito} \right) \\ & - \left( k_{alb,IL4}^+ C_{IL4} C_A - k_{alb,IL4}^- B_{IL4}^{alb} \right) \end{aligned} $$ $$ \frac{\partial B_{IL4}^{chito}}{\partial t} = k_{chito,IL4}^+ C_{IL4} (S_{chito} - B_{IL4}^{chito}) - k_{chito,IL4}^- B_{IL4}^{chito} $$ $$ \frac{\partial B_{IL4}^{alb}}{\partial t} = k_{alb,IL4}^+ C_{IL4} C_A - k_{alb,IL4}^- B_{IL4}^{alb} $$

Directly solving this complex system of three coupled PDEs and ODEs is computationally expensive and unnecessary. Considering that the molecular binding/dissociation reactions (nanosecond to millisecond scale) are much faster than the macroscopic diffusion process (hour scale), we employ the Quasi-Steady-State Assumption (QSSA) to simplify the model. This is a standard and effective method in receptor-ligand dynamics modeling (Lauffenburger & Linderman, 1993). We assume that the bound concentrations instantaneously reach equilibrium with respect to the free concentrations, i.e., \(\partial B / \partial t \approx 0\). This yields:

$$ B_{IL4}^{chito} \approx \frac{k_{chito,IL4}^+ S_{chito}}{k_{chito,IL4}^-} C_{IL4} := \alpha_{chito,IL4} C_{IL4} $$ $$ B_{IL4}^{alb} \approx \frac{k_{alb,IL4}^+ C_A}{k_{alb,IL4}^-} C_{IL4} := K_{A,IL4} C_A(x,t) C_{IL4} := \alpha_{alb,IL4}(x,t) C_{IL4} $$

The total concentration of IL-4 in the hydrogel is

$$C_{IL4}^{total} = C_{IL4} + B_{IL4}^{chito} + B_{IL4}^{alb} = (1 + \alpha_{chito,IL4} + \alpha_{alb,IL4}(x,t)) C_{IL4}$$

Applying the law of mass conservation to the total concentration, \(\partial C_{IL4}^{total} / \partial t = \text{Diffusion Term} + \text{Source/Sink Term}\), and rearranging, we obtain a modified, effective PDE that only involves the free concentration \(C_{IL4}\):

$$ \frac{\partial C_{IL4}}{\partial t} = \frac{1}{1+\alpha_{chito,IL4} + \alpha_{alb,IL4}(x,t)} \times \\ \left( D_{IL4,g} \frac{\partial^2 C_{IL4}}{\partial x^2} - (k_{IL4,g}^{deg} + k_{IL4,M} C_M) C_{IL4} + P_{IL4} \right) \quad $$

Similarly, for VEGF, we obtain a similar equation:

$$ \frac{\partial C_{VEGF}}{\partial t} = \frac{1}{1+\alpha_{chito,VEGF} + \alpha_{alb,VEGF}(x,t)} \times \\ \left( D_{VEGF,g} \frac{\partial^2 C_{VEGF}}{\partial x^2} - (k_{VEGF,g}^{deg} + k_{VEGF,M} C_M) C_{VEGF} + P_{VEGF} \right) \quad $$

In these two equations, the buffering term in the denominator, \(1 + \alpha_{chito} + \alpha_{alb}\), effectively retards the apparent diffusion, while the effective degradation rate in the parentheses, \(k^{deg} + k_M C_M\), quantifies the dynamic destructive effect of MMPs on drug stability.

Transport Model within the Tissue Domain (Ωt)

Once the therapeutic factors cross the interface into the tissue domain, \(0 \le x \le L\), they enter a pathological, heterogeneous environment dominated by granulation tissue. In the scenario of this study, the necrotic layer is approximately zero, and the tissue domain is divided into the granulation tissue region, \(0 \le x < x_g\), and the healthy tissue region, \(x_g \le x \le L\). Throughout the tissue domain, the factor concentrations are governed by the following general equations, where the diffusion and consumption parameters are piecewise constant according to the region:

$$ \frac{\partial C_{IL4}}{\partial t} = D_{IL4,t}(x) \frac{\partial^2 C_{IL4}}{\partial x^2} - k_{IL4,t}^{eff}(x) C_{IL4} \quad (5) $$ $$ \frac{\partial C_{VEGF}}{\partial t} = D_{VEGF,t}(x) \frac{\partial^2 C_{VEGF}}{\partial x^2} - k_{VEGF,t}^{eff}(x) C_{VEGF} \quad (6) $$

The diffusion coefficient in the necrotic zone, \(D_{i,n}\), is extremely low to reflect the dense tissue structure. In the granulation tissue region, which is the therapeutic target, the effective consumption rate, \(k_{i,gr}^{eff}\), is highest, lumping together all consumption processes, including cell receptor-mediated endocytosis. In the healthy tissue region, these parameters reflect baseline physiological levels.

Interface and Boundary Conditions

To obtain a unique solution for this PDE system, we impose well-defined initial and boundary conditions. At the initial time \(t=0\), the concentration of all species is zero. At the external boundaries of the system, \(x=-H\) and \(x=L\), a zero-flux condition, \(\partial C_i / \partial x = 0\), is applied. At the crucial gel-tissue interface, \(x=0\), the concentrations of MMPs and albumin are set to constant values, \(C_{M,source}\) and \(C_{A,source}\), to simulate a continuous source from the exudate. For the therapeutic factors, two conditions are simultaneously satisfied at the interface: concentration partitioning equilibrium, \(C_i(0^-,t) = \lambda_i C_i(0^+,t)\) (where \(\lambda_i\) is the partition coefficient), and mass flux continuity, \(-D_{i,g} \partial C_{i,g}/\partial x = -D_{i,t} \partial C_{i,t}/\partial x\). The criterion for efficacy is "near-interface target achievement": at the monitoring point \(x=0.02\,\text{mm}\), the first onset of action is considered to occur when \(C_{IL4}\ge1.5\,\text{nM}\) and \(C_{VEGF}\ge1.0\,\text{nM}\) are met.

Parameter Selection and Solution Method

To numerically solve the established PDE system, it is necessary to assign specific numerical values to all physicochemical parameters in the model. The selection of these parameters is key to the predictive power of the model. We have defined a "baseline scenario" for simulation based on extensive literature review, estimations from fundamental biophysical principles, and considerations for future measurable wet-lab experiments. This represents our best current estimate of the real diabetic chronic wound smicroenvironment. All parameters are in SI units.

The parameter values we selected, their physical meaning, units, and the basis for their selection are detailed in the table below:

Among these, parameters such as the diffusion coefficient (\(D\)) were estimated based on protein molecular weight using the Stokes-Einstein equation and were corrected for the tortuosity and density of different media (hydrogel, necrotic tissue, granulation tissue)^[4]. Reaction kinetics constants, such as the MMP degradation constant (\(k_{M}\)) and the tissue consumption rate (\(k_{eff}\)), were estimated based on published cytokine half-lives and enzyme kinetics studies^[5]. The yeast production rate (\(P\)) is considered a key "design parameter," with its value representing the production capacity we expect from the engineered yeast. It is worth noting that these parameters, especially the constants related to biological reactions, are core targets for subsequent parameter sensitivity analysis and experimental validation.

Since this model is a complex system composed of coupled, nonlinear PDEs, it has no analytical solution. Therefore, we used the Python language and the Finite Difference Method algorithm to perform numerical simulations of the system. Specifically, we discretized the continuous spatial domain into a uniform grid and used a second-order central difference scheme to approximate the second-order spatial derivatives (\(\partial^2 C / \partial x^2\)). for the time evolution, we used the explicit forward Euler method. This combination, known as the forward-time central-space (FTCS) scheme, transforms the complex PDE system into a set of discrete algebraic equations that can be solved iteratively over time steps on a computer. We set the total simulation time to 36 hours and the time step to 0.1 seconds for the numerical solution. The entire solving process was implemented using the NumPy library in Python.

Results

Spatiotemporal Distribution of IL-4 & VEGF Concentration

Under the setting of "necrotic layer thickness is approximately 0," the spatiotemporal heatmaps of IL-4 and VEGF show a rapidly established concentration gradient near the interface. The heatmap for IL-4 concentration is as follows:

The heatmap for VEGF concentration is as follows:

Time to Drug Onset near the Interface

Using \(x=0.02\,\text{mm}\) as the monitoring point (representing a depth of 20 micrometers into the wound surface), the time-course curves are shown below:

The results show: IL-4 first exceeds 1.5 nM at approximately 0.32 h, and VEGF first exceeds 1.0 nM at approximately 0.16 h.

Analysis & Discussion

Quantitative Prediction and Analysis of Onset Time

This computational simulation reveals the microscopic mechanism behind the rapid onset of action of the engineered yeast-hydrogel system. At the near-interface monitoring point (\(x=0.02\,\text{mm}\), representing a depth of 20 micrometers below the wound surface), VEGF reaches its therapeutic threshold (1.0 nM) within approximately 9.6 minutes (0.16 h), and IL-4 reaches its threshold (1.5 nM) within approximately 19.2 minutes (0.32 h). This finding has significant clinical implications: compared to traditional single-dose protein therapeutics, this "living" dressing system can establish effective dual-factor therapeutic concentrations near the wound surface within half an hour of application, providing a theoretical basis for the early initiation of anti-inflammatory and pro-angiogenic effects.

Notably, the onset time for VEGF is shorter than that for IL-4. This seemingly counterintuitive phenomenon can be explained by the following mechanisms:

1. 1. Dominant Role of Buffering Differences: Although VEGF has a larger molecular weight (~38 kDa) than IL-4 (~15 kDa), resulting in a lower intrinsic diffusion coefficient, IL-4 is subject to stronger binding buffering by chitosan in the gel (\(\alpha_{chito,IL4}=1.5\) vs. \(\alpha_{chito,VEGF}=0.5\)). This significantly reduces the effective diffusion coefficient of IL-4, \(D_{eff} = D/(1+\alpha_{total})\).
2. 2. Synergistic Impact of Degradation Rates: The baseline degradation rate of IL-4 (\(k_{deg}=0.684\,\text{h}^{-1}\), corresponding to a half-life of about 1 hour) is higher than that of VEGF (\(k_{deg}=0.2304\,\text{h}^{-1}\), corresponding to a half-life of about 3 hours), which further weakens its accumulation rate in the near-interface region.
3. 3. Partial Compensation by Production Rate: Although the production rate of IL-4 (\(P=36\,\text{nM/h}\)) is slightly higher than that of VEGF (\(P=28.8\,\text{nM/h}\)), this advantage is insufficient to counteract the combined adverse effects of buffering and degradation over the short time scale of half an hour.

Analysis of the Spatiotemporal Evolution of Drug Concentration

The process of drug concentration establishment can be clearly observed from the spatiotemporal heatmaps (Figure 1 and Figure 2):

• 1. Accumulation Phase within the Gel (0–10 h): In the initial phase of treatment, the continuously secreted IL-4 and VEGF first accumulate within the hydrogel layer (\(x<0\)). Due to the dual buffering effect of chitosan and albumin, the growth rate of the free concentration is effectively suppressed, and the concentration within the gel shows a near-linear increase in the first few hours, rather than instantaneous saturation.
• 2. Interface Penetration and Gradient Establishment (0.5–5 h): Once a concentration gradient is established within the gel, the drug molecules begin to diffuse through the interface (\(x=0\)) into the granulation tissue. The heatmaps show that observable concentrations appear near the interface (\(0 \le x \lesssim 0.1\,\text{mm}\)) within 0.5 hours, and a stable concentration gradient, decreasing from the interface into the tissue depths, forms over the subsequent hours.
• 3. Quasi-Steady State Achievement (10–36 h): After about 10 hours, the entire system enters a quasi-steady state, where continuous production in the gel, transport across the interface, and consumption within the tissue reach a dynamic equilibrium. At this point, the near-interface concentrations are maintained at a stable level above the threshold (IL-4 at ~2–3 nM, VEGF at ~1.5–2 nM), while concentrations at greater depths (\(x>0.5\,\text{mm}\)) are significantly lower due to tissue consumption and diffusion resistance.

Predictive Value and Limitations of the Model

This study has quantified the onset dynamics of the "living" dressing system through high-fidelity PDE simulations. Its predictive value lies in guiding future wet-lab animal experiments: without prior modeling, animal experiments might miss the optimal observation window by sampling too early (< 0.5 h) or too late (> 24 h).

However, the model also has certain limitations: as this simulation uses a linearized coupling method with "frozen coefficients"—treating MMP and albumin concentrations as known parameters at each time step—it may underestimate transient fluctuations in the presence of strong nonlinear feedback. Future work could incorporate inner iterations or a fully coupled implicit solver to improve accuracy. Additionally, we have neglected the two-dimensional healing dynamics and lateral diffusion from the wound edges towards the center, which might lead to an overestimation of the concentration in the central region. Furthermore, some parameters (such as \(K_A\), \(k_M\)) are based on literature estimates rather than actual measurements. Inter-patient variability in their true values could lead to a fluctuation of ±50% in the onset time.

Conclusion

Through linearized coupled reaction-diffusion PDE simulations, we have successfully predicted the rapid onset characteristic of the engineered yeast-chitosan hydrogel system in a scenario where the "interface is the granulation layer": VEGF and IL-4 reach their therapeutic thresholds near the wound surface in approximately 10 minutes and 20 minutes, respectively. The simulation further reveals that the onset time is primarily controlled by the "effective buffering" and "production rate" within the gel, while MMP degradation more significantly affects the long-term steady state and deep tissue exposure. These quantitative insights provide a solid theoretical foundation for subsequent formulation optimization (reducing buffering, increasing production rates, introducing MMP inhibitors), future animal experiment design (early sampling, near-interface detection), and clinical translation.

References

1. [1] Armstrong, D. G., Boulton, A. J. M., & Bus, S. A. (207). Diabetic foot ulcers and their recurrence. New England Journal of Medicine, 376(24), 2367–2375.
2. [2] Eming, S. A., Martin, P., & Tomic-Canic, M. (2014). Wound repair and regeneration: Mechanisms, signaling, and translation. Science Translational Medicine, 6(265), 265sr6.
3. [3] Wysocki, A. B., Staiano-Coico, L., & Grinnell, F. (1993). Wound fluid from chronic leg ulcers contains elevated levels of metalloproteinases bound to α2-macroglobulin. Journal of Investigative Dermatology, 101(1), 64–68.
4. [4] Saltzman, W. M. (2009). Drug delivery: Engineering principles for drug therapy. Oxford University Press.
5. [5] Rodriguez, J. P., González-Díez, M., Bustillo, A., González-Sarmiento, R., & Laso, F. J. (2006). Vascular endothelial growth factor is a substrate for matrix metalloproteinases-2 and -9. Biochemical and Biophysical Research Communications, 347(3), 798–804. https://doi.org/10.1016/j.bbrc.2006.06.143