Mathematical Prediction for Geometry‐Mediated Cell 3D In‐Growth on Bone Tissue Engineering Scaffolds
Xiang Gao, Zhijun Yu, Yu Yan, Jiamin Ding, Wangsiyuan Teng, Chenhe Zhou, Minjun Yao, Wenkan Zhang, Shenzhi Zhao, Fangqian Wang, Zhenxuan Shao, Jinfeng Zhou, Xiaoyong Wu, Chengcheng Yu, Liang Chen, Xiaoqiang Jin
TL;DR
This paper presents a mathematical model to predict how cell infiltration into bone tissue scaffolds is influenced by pore size and geometry.
Contribution
The study introduces a novel Porous-Fisher model that quantitatively predicts cell coverage rates and the impact of scaffold geometry on tissue growth.
Findings
Small pores promote horizontal bridging while large pores favor vertical migration of cells into scaffolds.
Convex geometries accelerate infiltration, while concave geometries allow spatiotemporal control of tissue growth.
Smaller pores in low diffusion environments may benefit elderly patients by delaying coverage.
Abstract
3D cell infiltration into porous scaffolds constitutes a fundamental prerequisite for bone tissue engineering. Though pore size and curvature are known to dictate this process, their mathematical coupling remains elusive. Herein, we identified a size‐dependent bone marrow‐derived mesenchymal stem cells 3D in‐growth pattern in which small pores promoted horizontal bridging, while large pores favored vertical cellular migration into the scaffold core. An analytical framework of Porous‐Fisher model was developed using a superposition approach tailored to boundary‐specific solutions. This approach not only enabled quantitative prediction of coverage rates through examination of grid dimensions and diffusion coefficients but also mathematically elucidated curvature and strategic geometric design. Furthermore, the prediction of cellular diffusion patterns on porous scaffolds was achieved…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5
FIGURE 6
FIGURE 7
FIGURE 8- —National Natural Science Foundation of China10.13039/501100001809
- —Key R&D Program of Zhejiang
- —Central Guidance on Local Science and Technology Development Fund of Zhejiang Province
- —5510 Project of the Second Affiliated Hospital of Zhejiang University
- —Natural Science Foundation of Zhejiang Province10.13039/501100004731
- —Zhejiang Administration of Traditional Chinese Medicine bureau and province co‐building major project
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsTopology Optimization in Engineering · Mathematical Biology Tumor Growth · Bone Tissue Engineering Materials
Introduction
1
Bone, a multifunctional tissue, plays a vital role in facilitating body movement and safeguarding internal organs. Bone defect, among the most prevalent traumatic injuries in humans, continues to pose significant challenges in contemporary clinical practice. Currently, autogenous bone graft remains the gold standard for clinical treatment, with annual statistics indicating that approximately 4 million people require bone transplantation or bone replacement surgery annually. Nonetheless, the donor sites scarcity, risk of secondary injury during harvesting, and challenge of achieving a precise anatomical fit are major hindrances to its widespread adoption [1, 2]. Therefore, these inherent constraints have necessitated more efficient strategies for bone regeneration.
The progress of tissue engineering has fostered the broad application of biomaterials in bone defect repair. Scaffolds with varied microarchitectures serve as temporary matrices for seeded cells, modulating their functions and guiding the regeneration of functional tissues and organs, thereby offering considerable potential in bone repair. Crucially, bone regeneration, a process encompassing angiogenesis, osteoconduction, and so forth, is largely driven by the orchestrated activity of multiple cell types [3]. Cell migration on the scaffolds, which underpins these biological processes [4, 5], is influenced by various biological and physical factors, including topological conditions. Despite recent advances, scaffold design remains at an early developmental stage [6]. While it is well‐established that cell can perceive such geometrical cues at subcellular scales and larger scales [7], but the underlying principles governing how scaffold geometry modulate cell behavior are still not fully understood. Thus, clarifying the relationship between geometry and cell migration is essential for designing optimized scaffolds adapted to diverse patient populations.
The geometry of porous scaffolds is defined by structural parameters such as pore shape, pore size, porosity, curvature, and so forth. Among these, pore size represents one of the most critical design variables. Current research suggests that effective bone formation requires pore sizes ranging from 300 to 1000 µm [8], with a minimum threshold of 100 µm, the precise criteria for tailoring pore size to specific pathological defects are currently lacking. Furthermore, local pore geometry undergoes dynamic changes during tissue development as a result of collective cell migration, which is itself guided by curvature‐dependent cues. The initial scaffold geometry can be described by pore size, whereas the evolving local geometry during cell migration is more accurately characterized by curvature. Notably, curvature exhibits an inverse relationship with pore size: larger pores correspond to lower average curvature values [9]. Despite extensive investigation into these phenomena, the mathematical relationship between geometry and migratory patterns of cells in varying states on the surface of scaffolds remains obscure. Thus, elucidating the interplay among pore size, curvature, and collective cell migration offers a promising direction for innovation in bone tissue engineering. A mathematical framework capable of predicting these interactions would significantly enhance scaffold design and application.
Herein, based on co‐culture experiments employing bone marrow mesenchymal stem cells (BMSCs) on porous scaffolds with varying pore size (400, 600, and 800 µm), the 3D in‐growth trend of BMSCs was systematically captured and analyzed (Figure 1). Notably, the migration trends were found to be closely associated with curvature. An analytical framework of the Porous‐Fisher model, derived via the superposition of particular solutions satisfying boundary conditions, not only facilitated quantitative prediction of cell coverage rates through examination of grid dimensions and diffusion coefficients but also provided mathematical insights into curvature and optimal geometric design. Furthermore, by modifying boundary conditions and diffusion coefficients, cellular diffusion patterns on porous scaffolds could be predicted. A deeper understanding of the mathematical and biological mechanisms governing cell ingrowth in relation to pore structure offers valuable theoretical guidance for the future clinical application of porous scaffolds.
The schematic illustrations of mathematical prediction and optimization for geometry‐mediated cell 3D in‐growth on cranial defect repairing scaffolds.
Experimental Section
2
Materials
2.1
Chitosan (CS) was purchased from Qingdao Haihui Bioengineering Co., Ltd (Biomedical grade, M_η_ = 9.1 × 10^5^ Da, D.D. = 89%). β‐TCP was obtained from Aladdin Reagent Inc. (Shanghai) and Sigma‐Aldrich, respectively. Other chemical reagents, including LiOH·H_2_O and urea, were all ordered from Sinopharm Chemical Reagent Co., Ltd. In addition, the live/dead cell staining kit was obtained from Biyuntian Biotechnologies (Shanghai, China). The cell counting kit‐8 (CCK‐8) kit, bicinchoninic acid (BCA) assay kit, and 4′,6‐diamidino‐2‐phenylindole (DAPI) were all purchased from Beyotime Biotechnology Co. (Shanghai, China). The primary antibodies of Collagen‐I, DyLight550‐labeled Donkey Anti‐Rabbit IgG (H+L) secondary antibodies, and Osterix (Ost) were bought from Saiye (Suzhou) Biotechnology Co., Ltd. (Suzhou, China), Boster Biotechnology Co. (Wuhan, China), and Affinity (China), respectively.
Preparation of Micropatterned Mold
2.2
To fabricate the micro‐pore porous scaffold, three micropatterned molds with distinct parameters were designed using AutoCAD commercial software (Autodesk, USA) [10]. As illustrated in Figure 2b, the size of molds was 12 mm × 12 mm and three types of micropillar arrays (400 µm × 400 µm micropillar array, 400 µm gap. 600 µm × 600 µm micropillar array, 400 µm gap. 800 µm × 800 µm micropillar array, 400 µm gap) were successfully processed by computer numerical control lathe (CNC, AD360 Ls, Sodick, Japan).
*(a) Schematic illustration of the fabrication of CS/β‐TCP scaffolds. (b) Design drawing of casting molds with different parameters. (c) Digital photos, pore models, and Micro‐CT images of CS/β‐TCP scaffolds (From top to bottom are CT400, CT600, and CT800). (d) Pore size parameters and (e) porosity of the scaffold material based on Micro‐CT results. (f) Relative BSA adsorption per unit area on hydrogel scaffold surface and (g) contact angle of the CS, CS/β‐TCP, CT400, CT600, and CT800 groups. (Statistical significance was set as ***p < 0.001, **p < 0.01, p < 0.05, and ns p ≥ 0.05).
Preparation of Micro‐Pore Porous Scaffold
2.3
The micro‐pore porous scaffolds were prepared via a two‐step method involving freeze‐thawing process and in situ gelation. First, 3 wt.% chitosan (CS) powder and 2 wt.% β‐tricalcium phosphate (β‐TCP) were dispersed in an 8.4 wt.% LiOH·H_2_O/6 wt.% urea solvent system through cyclic freeze‐thaw processing (−80°C freezing and 25°C thawing) [11]. The resulting CS/β‐TCP solution was then degassed and poured into the micropatterned mold. Three distinct CS/β‐TCP scaffold variants with different pore parameters were fabricated by gelation at 60°C for 4 h, followed by dialysis.
The CS/β‐TCP scaffolds with different pore parameters were named as CT400 (400 µm × 400 µm pore, 400 µm gap), CT600 (600 µm × 600 µm pore, 400 µm gap), and CT800 (800 µm × 800 µm pore, 400 µm gap) according to the pore size. Free of micropatterned mold, CS/β‐TCP hydrogels without pore were also prepared.
Characterization of Materials
2.4
The samples were lyophilized after frozen in −80°C and gold‐sprayed for conductance. The X‐Ray Diffraction (XRD, AXS D8 Advance X‐ray, Bruker, USA) spectra were recorded in the 2θ range of 5°–80° at a scan speed of 2°/min [12]. Image J was employed to analyze the CT‐reconstructed images of the scaffold for calculating pore size and porosity. The water contact angles of each material group were measured using a contact angle goniometer (OSA200‐T). To evaluate protein adsorption capacity, the samples (CT400, CT600, and CT800, 1cm×1cm×1 mm) were immersed in a 1% BSA solution and incubated at 37°C on a shaker for 8 h. After being rinsed 2–3 times with PBS, the samples were placed in a PBST solution and subjected to ultrasonication (100 W, 5–10 s per cycle, repeated ten times) at 4°C to elute the adsorbed BSA. The protein concentration in the supernatant was subsequently determined using a BCA assay kit (Beyotime, P0010S).
Biological Evaluation
2.5
Cytocompatibility and Proliferation Effect of Scaffold
2.5.1
The cytocompatibility and proliferation effect of scaffolds were evaluated through the CCK‐8 assay [13]. The human BMSCs (Cell Bank of Shanghai Institute of Biochemistry and Cell Biology) were selected as the model cell type due to their critical role in bone regeneration. To prepare the extracts, 100 mg CS/β‐TCP scaffold sample was incubated in 1 mL Dulbecco's modified Eagle's medium (DMEM), supplemented with 100 U/mL penicillin and 100 µg/mL streptomycin, for 24 h at 37°C [14]. The extracts were sterilized by filtration through a 0.22 µm membrane. BMSCs at a density of 5 × 10^3^ cells/well were cultured with 0.2 mL extracts in a 96‐well plate. After incubation at 37°C and 5% CO_2_ for 1, 3, and 7 days, the extracts were removed and CCK‐8 solution was added for additional 4 h’ incubation [15]. The absorbance of each well was measured by a micro‐plate reader at 450 nm.
Induction of BMSCs with Differential Migration Rates
2.5.2
To establish differential migratory capacities, BMSCs with different passages (P2, P6, and P10) were prepared by continuous passaging until P10 in DMEM high‐glucose medium supplemented with 10% FBS [16]. These cells were then seeded in 12‐well plates and cultured until full confluence. For the wound‐healing assay, the medium was replaced with serum‐free DMEM high‐glucose medium. A standardized scratch wound was subsequently created using a sterile 200 µL pipette tip. Cell migration dynamics were imaged over a 48 h period at 12 h intervals using an inverted phase‐contrast microscope (Olympus CKX53, Japan), with baseline wound dimensions recorded immediately post‐injury.
3D In‐Growth of BMSCs on Porous Scaffolds
2.5.3
The 3D cell ingrowth capacity into porous scaffolds was assessed using a co‐culture model [17]. Prior to cell seeding, all microporous scaffolds (radius = 6 mm, thickness = 1 mm) were sterilized by immersion in 75% ethanol (2 × 30 min), followed by thorough washing with phosphate‐buffered saline (PBS). The scaffolds were then transferred to an ultra‐low attachment 6‐well plate (Corning, Costar) and seeded with P2 BMSCs at a density of 2.5 × 10^5^ cells per well [18]. Cells were maintained under standard culture conditions (as described above), with medium replenishment every 48 h. After 1, 3, and 7 days of co‐culture, cell‐scaffold constructs were stained with Calcein AM/Propidium iodide (PI) and visualized using fluorescence microscopy (490/535 nm excitation/emission; STELLARIS5 confocal microscope, Leica, Germany) [19]. The experiment was repeated using senescent BMSCs to investigate the impact of cellular states on migratory behavior.
Cranial Defect Regeneration Studies In Vivo
2.5.4
All the Animal experiments were carried out under the approval of the Animal Experimentation Ethics Committee of the Second Affiliated Hospital of Zhejiang University School of Medicine (No. AIRB‐2023‐0784). Eight‐week‐old and eighteen‐month‐old male Sprague Dawley rats from the Academy of Medical Sciences of Zhejiang province were used. 3 days prior to surgery, the rats were transferred into the laboratory animal housing for acclimatization. Under general anesthesia, a 6.0 mm diameter cranial defect was produced using a trephine bur (external diameter 6.95 mm; internal diameter 6.0 mm) [20, 21, 22]. The sterilized cylindrical CS/β‐TCP hydrogel, CT400, CT600, and CT800, with a size of 6.0 mm in diameter and 1.0 mm in thickness (Each experimental group was divided into two subgroups: young and aged rats, with six rats per subgroup.) were implanted into the surgical sites, respectively. After a 4‐week postoperative period, euthanasia was performed by CO_2_ asphyxiation. The cranial specimens were collected and fixed using 4% paraformaldehyde for 48 h. Rats receiving no implant served as blank controls (defect‐only group).
Micro‐Computed Tomography
2.5.5
The condition of the newly formed bone at the cranial defect site 4 weeks after surgery, was scanned via Micro‐Computed tomography (Micro‐CT, Micro‐CT NV, Bruker, USA) [23]. Afterward, the SkyscanNRecon software was employed for 3D image reconstruction. All samples’ ROI was measured for bone volume (BV, mm^3^), total volume (TV, mm^3^), trabecular separation (Tb.Sp, mm), and the ratio of the segmented BV to the total volume of the region of interest (BV/TV, %).
Histological and Immunohistochemical Analysis
2.5.6
Following fixation in 4% paraformaldehyde and decalcification in 0.5 mol/L EDTA solution for 4 weeks, the specimens were paraffin‐embedded and sectioned at 5 µm thickness. Tissue sections were subjected to hematoxylin and eosin (H&E) staining for histological evaluation [24]. In addition, Masson and immunofluorescent staining were adopted to detect the detailed morphology of newly formed bone and expression of specific proteins using decalcified samples. Finally, the slices were photographed by microscope (BX53, Olympus, Japan) and fluorescence microscope.
Statistical Analysis
2.6
All the quantitative data were presented mean ± standard deviation (SD). Statistical analysis was carried out using One‐way ANOVA to determine the degree of significance. Statistical significance was set as ***p < 0.001, **p < 0.01 and *p < 0.05.
Results and Discussion
3
3D In‐Growth of BMSCs on Porous Scaffolds
3.1
Pore size is an important factor influencing bone in‐growth. Normally, a wide range of average pore sizes from 300 to 1000 µm is required to promote new bone formation [8]. To uncover the underlying mathematical relationship between cell 3D in‐growth and pore size of porous scaffolds, we engineered porous scaffolds with a pore size of 400, 600, and 800 µm via the casting method and specially‐made micropatterned molds (Figure 2a,b). The macroscopic and microscopic morphology were recorded in Figure 2c,d. The scaffolds exhibited intact structures with well‐defined geometries. The measured pore sizes of lyophilized CT400, CT600, and CT800 were about 420.18 ± 10.81, 620.93 ± 12.57, and 811.164 ± 16.95 µm, respectively. Further evaluation of hydrophilicity and protein adsorption behavior was conducted through contact angle measurements and protein adsorption assays. The results demonstrated comparable hydrophilicity and protein adsorption capacity among the three scaffolds (Figure 2f,g).
The 3D cell ingrowth capacity into porous scaffolds was evaluated using a co‐culture model with an ultra‐low attachment 6‐well plate. Fluorescence imaging of BMSCs after co‐culture with CT400, CT600, and CT800 for 1, 3, and 7 days showed that cell death was almost completely absent in these experiments (Figure 3b). Additionally, all three pore structures promoted transverse bridging tendency of BMSCs. The diffusion area increased from 23.6 ± 0.9% (CT400), 15.8 ± 0.96% (CT600), and 14.5 ± 1.1% (CT800) on day 1 to 49.0 ± 2.1%, 36.2 ± 1.34%, and 28.7 ± 3.65% on day 3, reaching 96.2 ± 1.9%, 50.6 ± 0.77%, and 38.2 ± 3.02% by day 7, respectively (Figure 3c). The average diffusion rate of BMSCs on CT400, CT600, and CT800 was calculated as 13.6 ± 0.5, 7.9 ± 0.7, and 5.0 ± 0.6%/day (Figure 3e). In parallel, as recorded by confocal microscopy (Figure 3b), BMSCs also could grow longitudinally into the pores of porous scaffolds. After 7 days of co‐culture, the penetration depth increased significantly from 100.0 ± 5.4 µm (CT400) to 156.3 ± 11.0 µm (CT600) and 236.9 ± 13.8 µm (CT800, Figure 3d). The corresponding infiltration rates were 14.3 ± 0.8, 22.3 ± 1.6, and 33.8 ± 2.0 µm/day (Figure 3e). The longitudinally infiltration trend was inversely proportional to the transverse bridging behavior. These observations suggested that pore size directly influenced BMSC spatial distribution. While smaller pore sizes promoted horizontal cell bridging, larger pores facilitated vertical infiltration into the scaffold depth, indicating the need for an optimized trade‐off [18]. Our findings suggested that CT600 achieved this balance most effectively, enabling integrated lateral and longitudinal cellular expansion. However, further investigation was required to clarify the mechanisms driving this behavior.
*(a) Schematic illustration of BMSCs 3D infiltration in CS/β‐TCP porous scaffolds. (b) Pore bridging experiments of BMSCs on CS/β‐TCP porous scaffolds for 1, 3, and 7 days, and the longitudinally in‐growth trend of BMSCs in CS/β‐TCP porous scaffolds for 7 days. (c) The plane coverage area proportion of porous scaffolds by BMSCs for 1, 3, and 7 days. (d) The longitudinal growth depth of BMSCs in porous scaffolds for 7 days. (e) The 3D in‐growth rate of BMSCs on porous scaffolds for 7 days. (Statistical significance was set as ***p < 0.001, **p < 0.01, p < 0.05, and ns p ≥ 0.05).
The fluorescence images of BMSCs revealed a collective cell migration pattern from the pore periphery to center, though further mechanistic investigation was required. To better explore the pore closure process, fluorescence intensity was color‐coded utilizing Python software on version 3.11.0 and Imaris software on version 9.7.2. As shown in Figure 4a, the migration front displayed higher fluorescence intensity, consistent with enhanced cell proliferation at the advancing edge and a preferential accumulation of BMSCs in this region [25]. Furthermore, the area inside the pore was extracted using the OpenCV package in Python 3.11.0, followed by 3D reconstruction using Matplotlib in Python 3.11.0 (Figure 4b,c). This analysis demonstrated three key findings: (1) Progressive expansion of BMSC‐infiltrated regions was observed, and the coverage rate of CT400 was highest in all groups. (2) Pore geometries exhibited gradual rounding, indicating preferential cellular colonization originating from angular regions. (3) Above 3D reconstruction results enabled clear correlation between the cell collective migration pattern with curvature, which was the “language of shape” [26]. However, curvature‐sensing mechanisms remained poorly characterized [27].
(a) The color‐coded fluorescence intensity of BMSCs collective migration and corresponding 3D reconstruction images. (b) Schematics of how to obtain projection and depth profile. (c) Front and left projection of the longitudinally in‐growth trend of BMSCs in CS/β‐TCP porous scaffolds for 7 days. Heatmap and depth profile corresponding to above confocal images.
In addition to top‐view observation of lateral cell ingrowth, side‐view projection was also employed to further analyze longitudinal cell penetration. Based on the frontal and left projection of the longitudinally in‐growth trend of BMSCs in CS/β‐TCP porous scaffolds for 7 days, cellular accumulation was predominantly localized in the upper pore regions, whereas deeper infiltration was observed in both CT600 and CT800 groups. These distribution patterns were further visualized through heatmap analysis using Python 3.11.0 and Imaris 9.7.2, highlighting preferential apical accumulation in CT400 scaffolds, where high‐intensity regions were concentrated in the pore throat. In contrast, BMSCs in the CT800 group primarily infiltrated longitudinally along the pore walls (Figure 3b), with overall cellular accumulation markedly lower than in CT400, as reflected in the more uniform distribution of signal intensity within the pores (Figure 4c). The CT600 group displayed intermediate characteristics between these two extremes. Depth profiling of longitudinal collective migration confirmed this trend. When combined with transverse bridging analysis, it was noted that complete transverse bridging only occurred in the CT400 group, while the achievement of transverse bridging might subsequently inhibit longitudinal infiltration, which might be attributed to the contact inhibition phenomenon [28], which was in accordance with our PCR verification (Figure S3). The aforementioned conclusion is further supported, indicating that the most balanced capability for promoting cellular infiltration was possessed by CT600.
Derivation of Cellular 3D In‐Growth Prediction Model
3.2
To investigate this phenomenon and establish the mathematical relationship between pore size, curvature, and BMSCs 3D in‐growth trends, a Porous‐Fisher model was employed in this study, with the analysis focused on a single pore. Although the Porous‐Fisher equation typically incorporated both diffusion and proliferation terms, the current investigation concentrated on diffusion dynamics within an individual porous structure. This simplification was justified as proliferative effects were considered negligible during BMSCs osteogenic differentiation when compared with diffusive processes [29]. The diffusion mechanism was governed by the standard diffusion equation:
where u(x, y, t) and α^2^ represented the cell density at point (x, y) at time t and diffusion coefficient, respectively. In this study, the following conditions were assumed for a square pore of side length L: (1) the cell density at the pore boundary was maintained at a constant value u_0_ (typically corresponding to the maximum density), and (2) the initial cell density within the pore interior (t = 0) was set to zero. These boundary and initial conditions were expressed mathematically as:
Through the variable transformation v(x, y, t) = u(x, y, t) − u 0, where the boundary value u_0_ was subtracted from u, zero boundary conditions were established for the new variable v(x, y, t), thereby simplifying the solution process. The solution was assumed to be separable and was expressed in the form:v(x, y, t) = X(x), Y(y), and T(t). When this expression was substituted into the diffusion equation and the variables were separated, a set of ordinary differential equations was obtained, each containing only a single independent variable. For the spatial component, since v was required to be zero at the pore boundaries, i.e., X(0) = 0, X(L) = 0 (similarly, Y(0) = 0, Y(L) = 0), The eigenfunctions satisfying these boundary conditions were found to be sine functions, with the wave numbers restricted to specific discrete values: sin(mπxL),sin(nπyL), Where m and n were positive integers. For a given spatial wave number pair (m, n), the temporal component satisfied: dTdt=−α2λmnT. This led to an exponentially decaying solution of the form: T(t)∝e−α2λmnt, where λ_mn_ was determined by the corresponding spatial wave numbers.
Utilizing the linearity property of the diffusion equation, a general solution was constructed by superimposing all valid particular solutions that satisfy the boundary conditions:
Owing to the symmetry of the initial and boundary conditions, only odd‐numbered terms for both m and n were required in the practical analysis, as given by:
In the early stages, before the cells had fully diffused throughout the system, BMSCs remained predominantly concentrated at the pore edges, where diffusion originated. As a result, higher‐order modes (corresponding to faster oscillations) dominated the initial dynamics. However, over time, these higher‐order modes decayed rapidly due to the exponential attenuation in the polynomial expression. After a sufficiently long period, only the lowest‐order mode (i.e., m = n = 1) retained a significant contribution, with its decay rate governed by:
Consequently, the primary characteristic time scale for the system to approach the boundary constant u 0 was:
When t ≫ τ_ m = 1, n = 1_, the cell density u(x, y, t) asymptotically approached u 0. Hence, if only a macroscopic estimation of the time required for the system to be effectively saturated was needed, retaining solely the contribution of the lowest‐order mode (m = n = 1) provided a reasonable approximation.
To further validate the theoretical analysis, the 2D diffusion equation was solved numerically using an explicit finite difference method. Additionally, finite element simulations were performed to model cell distribution at multiple time points. This approach enabled visualization of the diffusion process within the porous structure and verification of the analytical solution's accuracy.
The fundamental diffusion Equation (1) was expanded in a 2D Cartesian coordinate system as:
The spatial and temporal derivatives were discretized using the finite difference method. For the second‐order spatial derivatives, the central difference scheme was adopted:
The temporal derivative was discretized using the forward difference scheme:
By substituting these discretized forms into the original diffusion equation, the time‐stepping formula was derived in the following concise form:
For conciseness, the central point (i, j) was denoted as *u_p_ *, while its four adjacent points to the east, west, north and south as *u_e_ *, *u_w_ *, *u_n_ *, and *u_s_
- were similarly abbreviated. Using this notation, the governing equation was reformulated as follows:
To account for the complexity of actual cellular behavior, gradient terms were incorporated into the complete discretized formulation, yielding:
where the first term represented the contribution of the squared gradient and the second term represented the standard Laplacian operator. As the study primarily focused on standard diffusion processes, particular emphasis was placed on the Laplacian term in the numerical simulations.
To assess the accuracy of this model, 3D still frames of BMSC infiltration into scaffolds with varying pore architectures were simulated at 1‐, 3‐, and 7‐day intervals (Figure S5). The model‐generated outputs exhibited strong concordance with empirical observations of cellular migration patterns (Figure 3b). This close agreement between in silico predictions and experimental findings was observed between the simulated and experimental data, demonstrating that the model can be effectively applied to the simulation and prediction of 3D cellular ingrowth in porous scaffolds.
Prediction of Cell Migration Trend on Scaffolds Based on Model
3.3
The effect of grid size on diffusion dynamics was investigated numerically using square domains with side lengths of 400, 600, and 800 µm. A constant diffusion coefficient (α ^2^ = 0.25 µm^2^/time unit) was maintained throughout the investigations. The computational domains were discretized using a uniform spatial step size of 1 µm, resulting in grid configurations of 40 × 40, 60 × 60 and 80 × 80 cells, respectively. For all simulations, constant density initial boundary conditions (u = 1) were imposed on all four edges, while the initial density within the domain was set to 0. A fixed time step of Δt = 1 time unit was employed for all simulations. For each grid configuration, the characteristic time was defined as the moment when the 0.8‐density contour became fully filled, yielding values of approximately 1289, 2951, and 5292 time units for the 40 × 40, 60 × 60, and 80 × 80 grids, respectively (Figure 5a). Upon attainment of complete transverse bridging, the characteristic time threshold (1289 time units in this instance) was employed to model 3D cellular infiltration dynamics. As evidenced in Figure 5a, under uniform diffusivity (α ^2^ = 0.25 µm^2^ per time unit), the CT400 scaffold exhibited faster 2D bridging kinetics but showed the shallowest longitudinal infiltration depth—a consequence of contact inhibition. Thus, CT600 was determined to yield optimal 3D infiltration performance under these experimental conditions. The characteristic time demonstrated quadratic dependence on the linear dimension of the domain, in agreement with theoretical predictions for diffusion processes.
3D Cellular diffusion simulations: effects of boundary conditions and diffusion coefficients. (a) Boundary Condition Variation: cellular diffusion patterns on a 2D grid with different boundary dimensions (400, 600, and 800 µm) at a constant diffusion coefficient (α 2 = 0.25 µm2/time unit) and cellular diffusion patterns on a 3D grid at the characteristic time of 1289 time units. (b) Diffusion Coefficient Variation: cellular diffusion with different diffusion coefficients (α 2 = 0.125, 0.25, and 0.5 µm2/time unit) on a fixed boundary size of 600 µm and cellular diffusion patterns on a 3D grid at the characteristic time of 1476 time units.
To examine how the diffusion coefficient influences cellular migration patterns, a series of simulations on a fixed 60 × 60 grid (L = 600 µm) was performed. Three values of the diffusion coefficient (α ^2^ = 0.125, 0.25, and 0.5 µm^2^/time unit) were examined while maintaining consistent initial boundary conditions: constant density (u = 1) on all four edges and zero initial density within the domain. A fixed time step of Δt = 1 time unit was employed throughout. The simulations yielded characteristic times of approximately 5902, 2951, and 1476 time units for the respective diffusion coefficients (Figure 5b). This inverse relationship demonstrated that the characteristic time approximately halved when the diffusion coefficient doubled, consistent with fundamental diffusion theory, where the rate of density change was directly proportional to the diffusion coefficient. Besides, according to the cellular diffusion patterns on a 3D grid, under identical pore size (CT600), a marked enhancement in longitudinal wall‐adherent infiltration was observed in BMSCs under conditions of diminished diffusion coefficients (from 0.5 to 0.125 µm^2^/time unit). This approach enabled systematic examination of how diffusion speed influenced spatial distribution patterns of cellular density across the domain.
Effect of Diffusion Coefficients on 3D In‐Growth Patterns
3.4
Diffusion coefficients, recognized as a critical constant within the model, are determined by many factors, including cellular states and material properties. The diffusion coefficient is fundamentally defined as quantifying the cellular diffusion rate. To explore the influence of diffusion coefficients on cellular diffusion patterns, a scratch wound healing assay was performed on confluent cell monolayers using naive and senescent BMSCs. In our study, naive BMSCs exhibited a mean migration speed of 25.76 µm/h. In contrast, BMSCs at different stages of senescence (designated BMSCs(old) and BMSCs(older)) demonstrated significantly reduced migration velocities of 18.23 and 10.84 µm/h, respectively (Figure 6b). As the wound gap was completely closed by naive BMSCs within 48 h, migration rates were calculated at defined intervals (0–12, 12–24, and 24–36 h) to enable time‐resolved comparisons). Statistical analysis revealed that the migration rate of naive BMSCs was approximately 1.41 and 2.38 times higher than that of BMSCs(old) and BMSCs(older), corresponding to a diffusion coefficient ratio of roughly 1.99 and 5.64. Given the established relationship between cell migration velocity and diffusion coefficient in our model, these ratios could be proportionally applied (The anomalous diffusion exponent (α ^2^) was determined to be 0.25 for naive BMSCs, and thus, that of BMSCs(old) and BMSCs(older) were α ^2^ = 0.126 and 0.040 µm^2^/time unit). 3D Simulation results were visualized at multiple temporal stages for early‐stage observations with diffusion coefficients of α ^2^ = 0.126 and 0.040 µm^2^/time unit (Figure 6c,d). The simulations yielded characteristic times of approximately 2558 (400 µm, α ^2^ = 0.126), 5856 (600 µm, α ^2^ = 0.126), 10 500 (800 µm, α ^2^ = 0.126) and 8056 (400 µm, α ^2^ = 0.040), 18 444 (600 µm, α ^2^ = 0.040), 33 073 (800 µm, α ^2^ = 0.040) time units for the respective pore size, demonstrating that a delayed cellular coverage of porous bone repairing scaffolds was induced by lower diffusion coefficients. As shown in Figure 5a, the complete transverse bridging of naive BMSCs (α ^2^ = 0.25 µm^2^/time unit) was achieved at the characteristic time of 1289 time units, and 3D simulations were conducted at this time threshold. Comparative analysis revealed that at this defined time point, BMSCs with reduced diffusion coefficients (α ^2^ = 0.126 ‐ 0.040 µm^2^/time unit) failed to achieve complete planar coverage across CT400, CT600, and CT800 scaffolds, indicating the absence of contact inhibition. Thus, optimal performance in both transverse bridging and longitudinal infiltration was consistently observed in the CT400 group, and scaffold materials with smaller pore sizes might be more suitable for elderly bone defect patients. In clinical contexts, cellular infiltration capacity into porous scaffolds was significantly compromised by the aging process. Cellular populations derived from young and elderly patients were typically observed to exhibit distinct migration rates, which were matched with the above experimental scenario. These findings provided robust validation for the model's efficacy in predicting cellular infiltration patterns across scaffolds with varying pore geometries, while establishing its clinical utility for scaffold‐based therapeutic applications.
(a) Schematics of the influence factors for diffusion coefficient, and the study investigated how variations in the migration coefficient influence pore coverage patterns in scaffolds, with implications for clinical utility. (b) Scratch wound healing assay comparing naive and senescent BMSCs (BMSCs(old) and BMSCs(older)). 3D computational simulations (400, 600, and 800 µm domain) using diffusion coefficients α 2 = 0.126 µm2/time (BMSCs (old), (c)) and 0.040 µm2/time unit (BMSCs(older), (d)), demonstrating how experimental migration rates translate to diffusion model parameters.
Analysis of cellular diffusion in bone repair materials revealed that the propagation velocity of the diffusion interface was strongly correlated with local geometric characteristics. A distinct morphological transition was observed in the diffusion front, evolving “square to rounded rectangle to circular” morphology. To gain deeper insight into this phenomenon, the relationship between diffusion velocity and interface curvature was mathematically analyzed via the standard diffusion Equation (1):
An isoline u(x, y, t) = u 0 in the concentration field u(x, y, t) was defined as the interface of cellular diffusion. The normal velocity *v_n_
- of this interface (the rate at which the interface moves along its normal direction) could be expressed as:
After substituting the diffusion equation into the above expression:
In 2D space, the curvature κ of the interface was defined as:
Through mathematical analysis, the Laplacian operator could be decomposed into a form related to curvature:
where n=∇u|∇u| represented the unit normal vector of the interface. In the vicinity of the interface, gradient variations were found to occur predominantly along the normal direction, while changes in gradient magnitude at the interface itself remained relatively minor. Therefore, it could be approximated as:
Upon substitution of the aforementioned approximation into the normal velocity expression, the following relationship was obtained:
This key result indicated that the normal velocity of the interface was proportional to the local curvature, with the proportionality coefficient being the negative value of the diffusion coefficient.
For more general cases, where variations in gradient magnitude near the interface and additional influencing factors were considered, the normal velocity could be expressed as a curvature expansion:
where c0,c1,c2,… were coefficients determined by both the particular diffusion mechanism and local geometric characteristics. In many practical applications, sufficient accuracy was achieved through linear approximation, expressed as:
This theoretical framework elucidated the cellular behaviors observed in Figures 3 and 4 and mathematically clarified the underlying curvature‐sensing mechanism: the relationship between interface curvature and evolution velocity was deeply investigated through analysis of interfacial changes between consecutive time steps. For each angular correspondence, the Euclidean distance between matched points was calculated, representing the normal displacement of the interface. This displacement data was then correlated with the curvature distribution obtained from the preceding time step. Multiple regression analyses, including linear fitting, polynomial fitting, power‐law fitting, and exponential fitting were performed on the relationship between curvature and normal distance. The strength of curvature‐velocity correlation was quantitatively assessed through calculation of both Pearson correlation coefficients (measuring linear dependence) and Spearman rank correlation coefficients (evaluating monotonic relationships).
Role of Curvature on Cell Migration on Scaffolds
3.5
Numerical simulations were conducted on a 60 × 60 lattice (Lx = 60, Ly = 60) with an area fraction parameter of 0.25. The simulation ran for 3000 steps, with interface data collected at intervals of 60 steps. At each sampling point, the following measurements were stored for further analysis: Coordinates of all interface points and computed curvature values for each point. The evolution dynamics were analyzed by comparing consecutive time steps, where the Euclidean distance between corresponding points was calculated to determine normal velocity. Based on the simulation results, the diffusion interface evolution process was categorized into three distinct phases:
During the early stage (Figure 7a), the square pore boundary exhibited an approximately linear configuration with curvature values predominantly clustered near zero. Elevated curvature was observed only in limited sharp corner regions. Statistical analysis revealed a weak but significant curvature‐velocity correlation (Pearson r = 0.63, p < 0.001; Spearman r = 0.21, p = 0.04). Linear regression produced the relationship v n = 1.02 + 1.28κ, although the model exhibited poor goodness‐of‐fit (R ^2^ = 0.41). These findings suggested that interface dynamics during this stage were primarily dominated by initial conditions, surface potential energy distribution, or boundary diffusion factors rather than curvature‐dependent effects. While the gradient term contribution was confirmed, its influence appeared restricted primarily to high‐gradient regions proximal to boundaries.
(a) Schematics of the curvature‐dependent interface evolution at different simulation stages. Each panel (b–e) corresponds to a specific time step comparison (Step 1–61, 1201–1261, 1501–1561, and 2401–2461 respectively) and contains three related visualizations: Left: Boundary visualization showing the cellular diffusion interfaces at consecutive time steps (colored lines), illustrating the morphological evolution of the diffusion front. Middle: Angular distributions of curvature (blue) and normal displacement (red) along 0°–90° interfaces. Right: Correlation analysis between curvature (x‐axis) and normal displacement (y‐axis) with five mathematical fitting models (Linear, Quadratic, Cubic, Power, and Exponential).
During the middle stage (Figure 7b,c), curvature effects reached their peak here. As the evolution progressed, significant blunting of sharp corner regions was observed, with curvature emerging as the dominant kinetic parameter. A distinct positive correlation was established between curvature magnitude and normal velocity, demonstrating the characteristic “high curvature → high velocity” relationship. The correlation between curvature and normal velocity was remarkably strong (Pearson r > 0.98, p < 0.001; Spearman r > 0.99, p < 0.001). Linear fitting results (v n = 0.24 + 2.19κ for steps 1201–1261. v n = 0.29 + 2.16κ for steps 1501–1561) indicated close agreement with theoretically predicted values, thereby validating our analytical framework. During this phase, the contribution of the gradient term further diminishes, with diffusion kinetics being predominantly governed by the Laplacian term.
During the late stage (Figure 7d), exemplified by time step 2821, the interface geometry had considerably smoothed, closely approaching a circular configuration, which resulted in a markedly reduced variation in curvature across the interface (ranging from ≈0.004 to ≈0.145, with a mean of ≈0.122). Consequently, as curvature (κ) became nearly constant along the interface, the corresponding normal velocity (*v_n_ *) exhibited minimal fluctuation. Although a statistically significant positive correlation between curvature and normal velocity was still observed (Pearson r ≈ 0.44, p < 0.0001; Spearman r ≈ 0.50, p < 0.0001), its strength was substantially weaker than that in the middle stage. This weakening of this correlation was primarily attributed to the limited variation in the curvature itself. When the independent variable((κ)) exhibited a very a very narrow range, its capacity to statistically explain the (minor) variations in the dependent variable (*v_n_ *) decreased, even if the underlying physical relationship persisted. Thus, although the linear fitting yielded (*v_n_
- = 2.054 + 0.233κ), the explanatory value of this simple linear model was notably lower in this phase. At this stage, while more complex models (e.g., higher‐order polynomials) might have provided statistically better fits to the data, this likely reflected the model's ability to capture minor fluctuations or secondary physical effects (such as the relative influence of the gradient term on the small velocity variations) rather than indicating a fundamental change in the dominant kinetic mechanism.
Obviously, the bridging fronts of BMSCs progressively rounded prior to closure, demonstrating enhanced diffusion rates in convex corner regions and reduced rates in concave corners. As diffusion progressed, rapid smoothing of high‐curvature regions was observed, ultimately resulting in interface flattening with progressively diminishing curvature effects. Cellular diffusion rates could be selectively modulated through precise control of pore geometry, particularly curvature distribution, thereby enabling optimization of tissue regeneration processes. For practical implementation, regions requiring rapid cellular infiltration could be designed with predominantly convex corner structures, whereas concave geometries could be employed in areas where delayed growth was advantageous. These results demonstrated that the curvature‐velocity relationship exhibited dynamic evolution throughout the diffusion process, providing critical insights for predicting interfacial development across diverse material systems and important implications for porous bone repairing scaffolds design.
In Vivo Cranial Defect Regeneration
3.6
To assess the predictive validity of the pore‐size optimization derived from the above mathematical simulation modeling, we thus performed in vivo evaluation in a rat cranial critical‐size defect model. In our research, two holes of 6 mm diameter were drilled into the skull of rats, and scaffold samples with different pore sizes were placed in the holes (Figure 8a). Next, the regenerated new bones at 4 weeks were scanned via Micro‐CT analysis (Figure 8b). The blank group did not show distinct bone regeneration after 4 weeks, while the bone quality of the CS/β‐TCP group developed to a certain extent. However, the new bones in CS/β‐TCP group were mainly located on the edge of the hydrogels due to the unmatched degradation rate, which would make it hard for BMSCs to infiltrate in the hydrogels and cause insufficient release of Ca^2+^ and PO_4_ ^3−^ at the same time. Whereas, the new bone formation of porous scaffolds was better than that of CS/β‐TCP group. From the top view of Micro‐CT in the young rat group, the new bone growth into pore was detected, and the cranial defect area was almost completely covered by interconnected newly formed bone network after 4 weeks’ implantation in CT600 group. Though the new bone formation happened in the central region of cranial defect, the defect area was partially covered in CT400 and CT800. The best bone repairing effect appeared in CT600 group as expected. From the side view of Micro‐CT, the width of CS/β‐TCP group was smaller than that of the blank. There was still no new bone in the central defect area. Analogously, intermittent fresh bone tissue filled the defect in CT400, CT600, and CT800 groups. Furthermore, the filling degree of CT600 was highest. The filling types of CT400 and CT800 were point‐type and dash line‐type, respectively, which were in accordance with the projection and heatmap of cell experiments in Figure 5. In the aged rat group, although all experimental groups (CT400, CT600, and CT800) showed superior bone repair trends compared to the blank and β‐TCP/CS groups, the overall repair rate was significantly slower than that in the young group. Moreover, the repair trends among the three experimental groups in the aged cohort differed from those observed in the young group. Specifically, the CT400 group displayed the most favorable repair outcome in the aged rat cranial defect model, which aligns with our earlier mathematical predictions, thereby providing corroborative evidence for the predictive accuracy of our prediction model. The quantitative measurements of Micro‐CT, including BV, BV/TV, Tb.N, and Tb.Sp, which represented bone volume, bone volume fraction, the number of trabecular bones, and the trabecular bone spacing were also conducted (Figure S6). Compared to those of blank and CS/β‐TCP hydrogel group, the BV, BV/TV, and Tb.N of CT600 group in young rats significantly developed. Meanwhile, the Tb.Sp of CT600 was less than that of other groups. Furthermore, corresponding with the results of Micro‐CT, a few parts of incompletely degraded porous scaffolds remained in the newly formed bone tissue, which directly proving that BMSCs could infiltrate in the matrix of CT600 in vivo and osteogenesis of cranial defect proceeded from the surface of CT600 into pores. In comparison to young rats, aged rats exhibited significant reductions in BV, BV/TV, and Tb.N. Among all aged groups, the CT400 group demonstrated the highest levels of BV, BV/TV, and Tb.N, along with the lowest Tb.Sp. These results collectively indicated that CT400 offered the most effective cranial defect repair in aged rats. In addition, H&E, Masson, and immunofluorescent staining were further used to observe the detailed morphology of new bone and expression of specific proteins using decalcified samples. As depicted in H&E and Masson staining images (Figure 8b), new fibrous and bone tissue gradually formed to bridge the defect area in porous scaffold groups. Among these groups, there were voids and stripes in the regenerated tissue of CT400 and CT800 groups. Particularly, in CT600 group, almost all cranial defects were almost confluent, and the thickness of new bone was consistent with surrounding primary bone after 4 weeks, demonstrating that the remarkable bone tissue engineering regeneration ability of CT600 [30]. After quantitative analysis of collagen and mineral tissue area via Image J (Figure S7), more new bone and collagen formed in CS/β‐TCP hydrogel and CT600 group compared to the blank group. According to the immunohistochemical staining results of COL I and OST at 4 weeks (Figure 8b; Figure S6), the fluorescence area in both the CS/β‐TCP hydrogel and CT600 groups was significantly greater than that of the blank control, with CT600 exhibiting the most pronounced signal. The expression level of these osteogenic gene markers was in agreement with above H&E and Masson staining results. Collectively, these results corroborate the predictive capability of our simulation strategy for pore geometry and early‐stage cellular infiltration. Notably, CT600 scaffolds substantially accelerated cranial defect repair in rats by enhancing BMSC migration and osteogenic commitment.
(a) Schematic illustration of rat cranial defect model construction process. (b) 3D reconstructed Micro‐CT images (for the young and old group), H&E and MASSON staining images, immunofluorescence staining images of OST and COL I expression change (for the young group) of newly formed bone after implantation for 4 weeks.
Strengths and Limitations
4
In summary, this study advanced the field beyond current literature in several key aspects [7, 31, 32, 33]. First, a novel 3D mathematical model was developed, enabling more accurate prediction of complex cell infiltration dynamics than previously possible with conventional 2D frameworks. Second, a multi‐scale experimental validation system was established, incorporating primary replicative senescent cells, in vivo cranial defect repair, and aged animal models, thereby providing a physiologically relevant and translational evidence chain often lacking in prior studies. Finally, a geometry‐guided scaffold optimization strategy was proposed and experimentally validated through the design of age‐adapted pore architectures, offering both theoretical insights and practical guidance for transitioning bone regenerative materials from geometry‐driven design toward personalized medical applications. However, several limitations should be acknowledged. The current model remained simplified in terms of boundary conditions and cellular proliferation dynamics, which might not fully capture complex cell‐scaffold interactions. Additionally, the study relied on a constrained sample size and limited animal models, warranting further validation across broader pore size gradients, cell types, and temporal scales. While diffusion coefficients served as a primary mechanistic focus, more comprehensive biological pathways needed to be explored in future work. Although animal experiments exhibited trends consistent with computational predictions, definitive quantitative relationships remained to be established, particularly given the inherently multifactorial nature of bone regeneration that exceeded the current modeling framework.
Conclusion
5
By integrating mathematical simulations with experimental validation, this study elucidated the mechanistic role of pore curvature in dictating the spatiotemporal dynamics of BMSC infiltration. The Porous‐Fisher‐based framework developed here enabled systematic modulation of boundary conditions and diffusion coefficients, accurately replicating and predicting cellular coverage kinetics. Crucially, our model revealed a distinctive biphasic infiltration pattern, transitioning from curvature‐accelerated migration to contact‐inhibited stagnation, thereby establishing geometry as a governing variable in scaffold engineering. The modeling framework established a geometry‐guided strategy for scaffold optimization. From a structural perspective, convex topographies were shown to accelerate cellular infiltration in targeted regions, while strategically positioned concave domains modulated local growth kinetics, especially in contexts where delayed integration was desirable. Furthermore, simulations revealed that reduced diffusion coefficients delayed scaffold colonization, indicating that smaller pore‐size materials might better suit elderly bone defect patients. Ultimately, this predictive paradigm facilitated the rational design of patient‐specific scaffolds, thereby advancing the clinical utility of porous scaffolds.
Conflicts of Interest
The authors declare no conflicts of interest.
Supporting information
Supporting File: advs73608‐sup‐0001‐SuppMat.docx.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1B. Xu , P. Zheng , F. Gao , et al., “A Mineralized High Strength and Tough Hydrogel for Skull Bone Regeneration,” Advanced Functional Materials 27, no. 4 (2016): 1604327.
- 2B. M. Watson , T. N. Vo , A. M. Tatara , S. R. Shah , D. W. Scott , P. S. Engel , and A. G. Mikos , “Biodegradable, Phosphate‐containing, Dual‐gelling Macromers for Cellular Delivery in Bone Tissue Engineering,” Biomaterials 67 (2015): 286–296, 10.1016/j.biomaterials.2015.07.016.26232878 PMC 4550499 · doi ↗ · pubmed ↗
- 3Y. Zhang , P. Wang , J. Jin , et al., “In Silico and in Vivo Studies of the Effect of Surface Curvature on the Osteoconduction of Porous Scaffolds,” Biotechnology and Bioengineering 119, no. 2 (2022): 591–604, 10.1002/bit.27976.34723387 · doi ↗ · pubmed ↗
- 4C. Leclech , D. Gonzalez‐Rodriguez , A. Villedieu , T. Lok , A. M. Déplanche , and A. I. Barakat , “Topography‐Induced Large‐Scale Antiparallel Collective Migration in Vascular Endothelium,” Nature Communications 13, no. 1 (2022): 2797, 10.1038/s 41467-022-30488-0.PMC 912015835589751 · doi ↗ · pubmed ↗
- 5R. Mayor and S. Etienne‐Manneville , “The Front and Rear of Collective Cell Migration,” Nature Reviews Molecular Cell Biology 17, no. 2 (2016): 97–109, 10.1038/nrm.2015.14.26726037 · doi ↗ · pubmed ↗
- 6H. Wu , X. Wang , G. Wang , et al., “Advancing Scaffold‐Assisted Modality for In Situ Osteochondral Regeneration: a Shift from Biodegradable to Bioadaptable,” Advanced Materials 36, no. 47 (2024): 2407040, 10.1002/adma.202407040.39104283 · doi ↗ · pubmed ↗
- 7S. J. P. Callens , D. Fan , I. A. J. van Hengel , et al., “Emergent Collective Organization of Bone Cells in Complex Curvature Fields,” Nature Communications 14, no. 1 (2023): 855, 10.1038/s 41467-023-36436-w.PMC 998448036869036 · doi ↗ · pubmed ↗
- 8M. N. Collins , G. Ren , K. Young , S. Pina , R. L. Reis , and J. M. Oliveira , “Scaffold Fabrication Technologies and Structure/Function Properties in Bone Tissue Engineering,” Advanced Functional Materials 31, no. 21 (2021): 2010609, 10.1002/adfm.202010609. · doi ↗