Abstract

Background: Adult orthodontic tooth movement occurs within a mechanically constrained alveolar environment, where bone remodelling capacity is limited and biological responses are spatially heterogeneous. A previous finite element study (Part I) characterised the mechanical environment of orthodontic tipping in adult alveolar cortical bone as a system of structured spatial gradients rather than uniform stress or strain distributions. While that description provided a mechanically rigorous framework, it remained inherently static and did not address how repeated orthodontic loading may interact with the adaptive capacity of mature alveolar bone over time.

Objectives: To extend the gradient-based mechanical framework established in Part I into a time-dependent, iterative mechanobiological formulation; to implement an automated Python post-processing pipeline for spatial gradient computation and harmonic mesh deformation; and to validate the computational stability and geometric robustness of this closed-loop FEA-Python architecture as a prerequisite for patient-specific software integration (Part III).

Materials and Methods: The finite element geometry described in Part I was reused as the baseline configuration. A gradient-based stimulus function Φ was defined as a weighted combination of hydrostatic stress, its spatial gradient and the difference between principal stresses. An iterative pipeline integrating ANSYS Workbench finite element analysis with external Python post-processing was implemented. Nodal mechanical variables were exported and indexed via a k-dimensional tree (KD-Tree) with a 1.5 mm search radius; spatial gradients were computed via inverse-distance-weighted least-squares regression. Geometric adaptation was resolved across 30 iterative computational remodelling cycles corresponding to a single aligner activation stage prescribing a total coronal displacement of 0.25 mm. Each iterative cycle represented a mechanical equilibrium step with incremental density updating (Δρ), rather than a sequential clinical aligner stage. A power-law density-modulus relationship and a linear stimulus-driven density update were used. Loading conditions were kept constant throughout iterations to isolate the effect of evolving material properties on the mechanical response.

Results: The iterative pipeline produced stable and convergent evolution of bone density and elastic modulus. Density variation was spatially organised, with localised increases in compression-dominated regions and relative reductions in tension-dominated regions, separated by smooth transitions rather than sharp boundaries. The Python post-processing module achieved 100% numerical stability across all benchmark cycles. In a control sample of 500 cortical elements under vestibular tipping, hydrostatic stress ranged from −0.0767 MPa (maximum focal compression at the cervical margin) to +0.1007 MPa (maximum tension on the opposing wall), consistent with the low-force regime reported in Part I. The mean spatial gradient magnitude was 0.0061 MPa/mm, of the same order as lacunar pressure differentials implicated in osteocyte mechanosensing. The harmonic deformation algorithm assimilated gradient-driven displacement vectors without mesh distortion, folding or inverted elements, enabling smooth anatomical morphing of the gingival envelope across all 30 increments.

Conclusion: Adult alveolar cortical bone adaptation under sustained orthodontic loading appears to be governed by the spatial organisation of the mechanical field rather than by absolute stress or strain magnitudes. The proposed iterative framework extends the static gradient-based interpretation of Wolff’s law established in Part I into a temporal mechanobiological formulation, while demonstrating that a closed-loop FEA-Python pipeline is computationally robust and geometrically stable. The successful integration of spatial gradient computation, harmonic mesh deformation and iterative material updating provides a mechanically consistent foundation for the subsequent translation of this framework into orthodontic treatment-planning software (Part III).

Keywords: Orthodontic Biomechanics; Wolff’s Law; Iterative Finite Element Analysis; Bone Remodelling; Mechanical Gradients; Mechanobiology; Adult Orthodontics; Hydrostatic Stress; Computational Pipeline; Harmonic Mesh Deformation

Introduction

Orthodontic tooth movement in adult patients is a biomechanical process that takes place in a mechanically constrained alveolar environment, where bone remodelling capacity is reduced and biological responses are highly context-dependent [1-4]. Within this setting, mechanical stimuli generated by orthodontic forces do not translate into uniform or predictable patterns of adaptation. Instead, they define spatially heterogeneous mechanical fields whose biological relevance remains conditional, localised and modulated by patient-specific factors [2-6]. This non-deterministic behaviour challenges classical mechanostat-based interpretations, which rely on threshold-driven relationships between strain magnitude and remodelling response and assume a relatively direct coupling between mechanical stimulus and biological outcome [1,7].

In a previous finite element study, the mechanical environment associated with orthodontic tipping in adult alveolar cortical bone was characterised as a system of structured spatial gradients rather than uniform stress or strain distributions [8]. Both coronal and radicular tipping configurations generated heterogeneous fields of strain, hydrostatic stress and principal stresses, with spatial gradients across cortical thickness and along the root surface defining mechanically distinct domains. Importantly, global strain magnitudes remained within low-to-moderate ranges consistent with an adult homeostatic context, while local gradients and stress transitions delineated zones of potential biological relevance without implying deterministic remodelling outcomes [8]. This field-based description provided a mechanically rigorous framework for understanding orthodontic loading; however, it remained inherently static and did not address how repeated or sustained orthodontic stimuli may interact over time with the adaptive capacity of mature alveolar bone. Moreover, the constitutive modelling of the periodontal ligament as a linearly elastic material in that initial framework did not capture the time-dependent stress relaxation that characterises collagenous soft-tissue response under sustained mechanical loading. The present study extends the static gradient-based description by incorporating a viscoelastic periodontal ligament formulation (generalised Maxwell-type: αiG = 0.925, αiK = 0.890, τi = 0.15 s), thereby allowing the mechanical field to redistribute across the alveolar cortex as the ligament relaxes over the characteristic time scale of orthodontic force application. This temporal load-transfer mechanism is particularly relevant in adult orthodontics, where sustained aligner-mediated displacement generates a quasi-static bone stimulus after initial PDL relaxation and where the spatial organisation of the resulting cortical stress gradients may dominate the biological response.

Bone adaptation under sustained mechanical loading has been the subject of extensive computational investigation. Iterative finite element approaches have been developed to simulate long-term orthodontic tooth movement by alternating mechanical computation with bone remodelling steps, in which tissue properties or geometry are updated in response to the calculated mechanical state [2,9-12]. Most of these formulations rely on strain energy density or strain-magnitude criteria as the driving variable, assuming a homeostatic equilibrium against which deviations trigger anabolic or catabolic responses [1,9,12]. Although such approaches have provided valuable insight into long-term tooth displacement, they typically assume a relatively direct stimulus-response coupling and do not incorporate the spatial organisation of the mechanical field as a primary descriptor of the mechanobiological environment.

Recent advances in bone mechanobiology indicate that osteocyte responses depend not only on the magnitude of mechanical loading but also on its spatial organisation and on local gradients of deformation and pressure [5,13-17]. Hydrostatic stress gradients in particular have been associated with interstitial fluid flow modulation, lacunar pressure build-up and osteocyte mechanosensing within the lacuno-canalicular network [5,13,16,17]. Mechanosensitive ion channels such as PIEZO1 have been identified in periodontal ligament cells and bone cells as plausible mediators of such mechanotransduction processes during orthodontic tooth movement, although their precise role in adult alveolar bone remodelling remains under investigation [3,4,7,18,19]. These findings — together with earlier in vivo evidence that strain gradients correlate with sites of periosteal bone formation, and with continuum models incorporating strain-gradient descriptors in cortical bone adaptation — support the hypothesis that, under adult orthodontic conditions, the spatial organisation of the mechanical field may play a more relevant role than absolute stimulus magnitude in defining regions of potential adaptation [2,7,8,13,20,29].

The aim of the present study was to extend the static gradient-based mechanical framework of Part I into a time-dependent, iterative mechanobiological formulation, in which a closed-loop FEA-Python pipeline drives sequential updates of bone density and elastic modulus through a spatially organised mechanical stimulus. By explicitly separating mechanical computation from biological interpretation, this framework is intended to validate the computational stability and geometric robustness required for its subsequent integration into patient-specific orthodontic treatment-planning software (Part III).

Materials and Methods

A time-dependent computational framework was developed to investigate the evolution of adult alveolar cortical bone properties under sustained orthodontic loading. The model integrates finite element analysis with an external Python-based iterative algorithm, enabling sequential updates of bone density and corresponding material properties. This formulation extends the static mechanical characterisation described in Part I into a dynamic context in which the interaction between mechanical stimuli and structural properties is evaluated across multiple simulation cycles [8]. The objective is not to reproduce biological remodelling processes in detail, but to explore how spatially organised mechanical fields may influence structural evolution within a simplified mechanobiological context [2,9-11].

The finite element model was derived from the geometry previously described in Part I and was used as the baseline configuration for all simulations [8]. The model represents a segment of the maxillary alveolar bone including the tooth, the periodontal ligament and the surrounding cortical and trabecular bone compartments. The periodontal ligament was modelled as a uniform layer of 0.2 mm thickness with a viscoelastic constitutive formulation to capture time-dependent stress relaxation and load redistribution under sustained displacement-controlled loading. The generalised Maxwell-type viscoelastic behaviour was parameterised by a relative shear modulus αiG = 0.925 (corresponding to 92.5% shear relaxation), a relative bulk modulus αiK = 0.890 (corresponding to 89% volumetric relaxation) and a characteristic relaxation time τi = 0.15 s. These values were selected to approximate the non-linear, strain-rate-sensitive response of collagenous periodontal fibres under orthodontic force levels in the low-force regime (2-3 N), consistent with the compliant mechanical environment characteristic of aligner-mediated tooth movement [2,11]. It should be emphasised, however, that the relaxation time τi = 0.15 s adopted here represents a numerical simplification chosen to ensure computational stability of the viscoelastic solver within a single load step, rather than a physiological clinical time scale; the biologically relevant relaxation times of the periodontal ligament under sustained orthodontic loading extend over minutes to hours and are here absorbed by the quasi-static equilibrium assumption enforced at each iteration cycle. The alveolar cortical bone thickness was set to 1.7 mm based on direct CBCT-derived measurements, as detailed in Part I [8]. A three-dimensional mesh composed of approximately 65,000 tetrahedral elements and 120,000 nodes was generated. The alveolar cortical and trabecular bone compartments were modelled as homogeneous, isotropic and linearly elastic, in line with commonly adopted assumptions in orthodontic finite element studies [2,8,11]. The periodontal ligament was assigned the viscoelastic formulation described above to represent time-dependent shear and volumetric relaxation, while the tooth (enamel-dentin complex) was treated as rigid relative to the surrounding soft tissues. All other structural components (aligner, attachment, gingival envelope) were modelled as linearly elastic with properties consistent with the static baseline configuration reported in Part I [8]. Although these simplifications do not capture the anisotropic and viscoelastic behaviour of biological tissues, they provide a controlled environment for analysing mechanical-structural interactions over successive iterations [1,2]. Contact interfaces between the tooth, periodontal ligament and bone were defined as bonded, while the contact between the aligner and the tooth surface was modelled as frictionless, allowing relative sliding while transmitting normal contact forces [8]. A fixed boundary condition was applied at the basal region of the model to prevent rigid body motion. Orthodontic loading was applied using a displacement-controlled approach consistent with aligner-mediated mechanics, prescribing a coronal displacement of 0.25 mm at the aligner-tooth interface. The resulting force system was not imposed directly but emerged from the contact interaction, generating a reaction force in the low-force regime (approximately 2-3 N) consistent with previously reported aligner-based orthodontic systems [8]. The loading configuration was kept constant throughout all iterations in order to isolate the effect of evolving material properties on the mechanical response of the system.

Following each finite element simulation, nodal mechanical variables were exported for external post-processing. The variables considered in this study included equivalent elastic strain (εeq), hydrostatic stress (σh) and the maximum and minimum principal stresses (S1 and S3), defined consistently with the descriptors introduced in Part I [8]. In addition to these variables, the spatial gradient of hydrostatic stress |∇σh| was computed numerically from the nodal stress distribution. This gradient term, defined as the magnitude of the spatial variation of σh, was incorporated as a key descriptor of mechanical heterogeneity within the alveolar cortical bone, in line with experimental and theoretical evidence indicating that osteocytes are sensitive to spatial differences in mechanical state rather than to homogeneous loading alone [5,13,14,16,17]. The extraction and processing of these variables were performed externally using Python, allowing flexibility in the definition of the mechanobiological stimulus and avoiding the incorporation of biological assumptions within the finite element solver [2].

A gradient-based stimulus function Φ was defined to describe the mechanical environment driving structural adaptation:

Φ = w₁·σ_h + w₂·|∇σ_h| + w₃·(S₁ − S₃)

where w₁, w₂ and w₃ are weighting coefficients that allow the stimulus function to incorporate volumetric stress states, spatial heterogeneity and directional asymmetry of loading. Unlike classical formulations based solely on strain magnitude or strain energy density, this formulation emphasises spatial transitions within the mechanical field rather than absolute values alone, in agreement with continuum bone remodelling models that explicitly distinguish hydrostatic and deviatoric components of deformation [7,9,10,12,13]. Bone density evolution was modelled using an incremental update rule applied at each iteration step:

Δρ = k · Φ · Δt

where k is a remodelling rate constant and Δt is the time increment between iterations. This formulation provides a simplified representation of mechanobiological adaptation, in which mechanical stimuli influence structural properties over time [9,10]. The update rule does not represent a biological law but serves as a computational mechanism to explore the coupling between mechanical fields and material structure. The updated density values were then used to modify the elastic modulus of the bone according to the standard power-law relationship widely adopted in computational bone mechanics [9,12,22].

E = E₀ · (ρ / ρ₀)ⁿ

where E₀ is the reference elastic modulus, ρ₀ is the initial density and n is an empirical exponent. This relationship introduces a feedback mechanism in which local changes in density alter the stiffness of the material, thereby influencing the distribution of mechanical stresses in subsequent simulations [9,10,21,22]. The computational framework followed a closed-loop iterative process consisting of: finite element simulation to compute the mechanical field, export of nodal mechanical variables, calculation of the stimulus function Φ in Python, update of bone density ρ, update of elastic modulus E, reassignment of material properties in the finite element model and execution of a new simulation cycle. This process was repeated for a predefined number of iterations, allowing evaluation of the temporal evolution of the system. Representative parameter values were selected to ensure numerical stability and physiologically plausible behaviour, with a remodelling constant k = 2 × 10⁻³ mm³/(MPa·day), an exponent n = 2.5 and a time step Δt = 1 day. These values were not calibrated against experimental data but were chosen to maintain a stable and controlled evolution of the model, consistent with parameter ranges previously reported for cortical bone adaptation studies [9,21,22]. Results were analysed by evaluating the spatial distribution of bone density and elastic modulus across iterations. Density maps were generated at selected time points to assess the evolution of structural patterns and qualitative spatial comparisons were performed between regions of high mechanical gradients and regions exhibiting significant density changes [20,29]. The objective of this analysis was to identify spatial correlations between mechanical heterogeneity and structural adaptation, without introducing biological thresholds or predictive criteria. The closed-loop iterative architecture integrating ANSYS Workbench finite element analysis with external Python post-processing is illustrated in Fig. 1.

Figure 1: Closed-loop iterative FEA-Python pipeline architecture. The cycle alternates between ANSYS Workbench (mechanical computation) and Python post-processing (gradient indexing, stimulus evaluation, density update). Nodal variables (εeq, σh, S1, S3) are exported and indexed via KD-Tree (search radius 1.5 mm); spatial gradients ∇σh are computed via inverse-distance-weighted regression. Density updates follow Δρ = k·Φ·Δt and modulus updates follow E = E0·(ρ/ρ0)ⁿ. Pipeline architecture adapted from iterative FEA-remodelling schemes reported in [9,12,22]. Parameters: k = 2×10−3 mm³/(MPa·day), n = 2.5, Δt = 1 day, w1 = 0.3, w2 = 0.5, w3 = 0.2, ρ0 = 1.8 g/cm³, E0 = 13.7 GPa.

To enable the integration of this iterative framework into orthodontic treatment-planning software (Part III), a dedicated Python post-processing module was developed for spatial gradient computation, nodal indexing and harmonic mesh deformation. The technical implementation of these geometric algorithms is detailed below.

Spatial Gradient Computation and Geometric Post-Processing

To translate the three-dimensional mechanical fields computed by the finite element solver into a format suitable for iterative visualisation and subsequent integration into orthodontic treatment-planning software (Part III), an automated Python post-processing module was implemented. Following each simulation cycle, nodal hydrostatic stress values (σh) were exported from the cortical bone domain and indexed in a spatial data structure. A k-dimensional tree (KD-Tree) was constructed over the centroid cloud of the cortical mesh elements to enable efficient neighbourhood queries. A critical search radius of 1.5 mm was selected to capture the local lacuno-canalicular interaction domain while preserving computational tractability, consistent with the spatial scales reported for osteocyte mechanosensing and lacunar pressure propagation [5,16,17].

For each volumetric element i, the spatial gradient of hydrostatic stress (∇σ_h) was approximated via an inverse-distance-weighted least-squares regression over the local neighbourhood:

∇σ_h⁽ⁱ⁾ ≈ ( Σ_j w_j · (σ_h^(j) − σ_h⁽ⁱ⁾) · d_ij ) / ( Σ_j w_j )

where d_ij represents the inter-centroid distance vector and w_j = ||d_ij||⁻² acts as a geometric penalty factor. This formulation yields a three-dimensional vector field representing the direction and magnitude of pressure flux across the alveolar cortex, consistent with the gradient-based mechanical descriptors introduced in Part I and with continuum frameworks in which hydrostatic-stress gradients are treated as a primary remodelling descriptor [8,13,20,29]. For geometric adaptation across the 30 aligner activation increments (0.25 mm per stage), a harmonic mesh-deformation algorithm was implemented. The displacement field was formulated as a Dirichlet boundary-value problem governed by the discrete Laplace-Beltrami operator (ΔU = 0). The Laplacian matrix L was assembled using cotangent weights to preserve anatomical smoothness and prevent mesh folding or geometric singularities:

L_ij = (1 / 2A_i) (cot α_ij + cot β_ij)

The resulting linear system L U = B was solved with the mechanical gradient vectors incorporated as deformation drivers at the mobile boundary nodes, yielding the updated three-dimensional vertex coordinates for each soft-tissue and cortical mesh. This geometric update step is performed externally to the finite element solver, maintaining the explicit separation between mechanical computation and biological or geometric interpretation advocated in the present framework [2,13,28,29].

For qualitative visualisation purposes within the post-processing interface, the hydrostatic stress field was partitioned into three descriptive bands: compression-dominant regions (σh < −0.02 MPa), homeostatic-equilibrium regions (−0.02 MPa ≤ σh ≤ 0.05 MPa) and tension-dominant regions (σh > 0.05 MPa). These bands are used for colour-mapping and clinical orientation only; they do not constitute biological remodelling thresholds and are not incorporated into the quantitative density update rule (Δρ = k·Φ·Δt).

The algorithmic implementation of the localised spatial neighbour queries was written in Python 3.8 using the SciPy library. The source code and processed datasets supporting the findings of this study are available from the corresponding author upon reasonable request.

Results

The implementation of the iterative finite element-Python pipeline resulted in a stable and reproducible evolution of bone density and elastic modulus across successive simulation cycles. All simulations were performed under constant boundary conditions, allowing isolation of the effects of material property updates on the mechanical field. The iterative process produced progressively modified mechanical environments, reflecting the coupling between local stiffness and load distribution. Progressive changes in bone density were observed across iterations. During early iterations, density changes were spatially localised and closely followed the regions of high mechanical gradient identified in the initial configuration described in Part I [8]. As iterations progressed, density redistribution became more structured, with coherent regions of increased and decreased density emerging across cortical thickness and along the root surface. In late iterations, the system approached a quasi-stationary configuration, with the magnitude of density updates between successive cycles decreasing progressively. Density variations remained within moderate ranges throughout the simulation, without abrupt or unstable transitions. The system did not diverge or exhibit oscillatory instabilities, in line with previous reports on the convergence behaviour of similar density-based remodelling algorithms [9,11,21,22]. The temporal evolution of bone density across successive remodelling cycles, from the homogeneous baseline to the quasi-stationary configuration, is shown in Fig. 2.

Figure 2: Temporal evolution of bone density (ρ) across iterative remodelling cycles. (A) Iteration 0: homogeneous baseline (ρ0 = 1.8 g/cm³). (B) Iteration 10: early localised changes in compression-dominated regions (green, cervical margin) and tension-dominated domains (purple, opposing wall). (C) Iteration 20: structured redistribution with smooth transitional zones. (D) Iteration 30: quasi-stationary configuration. Tooth root (grey) and PDL (blue dashed, 0.2 mm) are excluded from the remodelling map. Scale: mm. Colormap: PiYG.

Density evolution was not uniform. Instead, the model generated spatially organised remodelling patterns, characterised by localised increases in density in regions associated with compression-dominated stress states (S3 < 0), relative reductions in density in regions associated with tension-dominated domains (S1 > 0) and smooth transitions between these regions rather than sharp boundaries. These patterns were consistent with the spatial organisation of mechanical gradients rather than with absolute stress or strain magnitudes, in agreement with the heterogeneous and non-uniform behaviour of adult alveolar cortical bone reported under orthodontic loading conditions [2,5,6,8].

A clear spatial correspondence was observed between regions of high hydrostatic stress gradients (|∇σ_h|) and regions exhibiting the greatest density variation (Δρ). Specifically, zones with elevated gradient magnitude showed higher rates of density change, whereas mechanically homogeneous regions exhibited minimal or negligible density updates. This spatial relationship remained consistent throughout the iterative process. From a mechanobiological perspective, this is compatible with current evidence indicating that osteocyte responses are modulated by interstitial fluid flow and local pressure differentials rather than by homogeneous strain alone [5,13,16,17]. Changes in bone density resulted in corresponding updates of the elastic modulus distribution, with regions of increased density exhibiting local stiffening and regions of decreased density showing relative compliance. These changes altered the load transfer within the alveolar bone, producing secondary modifications in the mechanical field at subsequent iterations. The iterative framework therefore revealed a bidirectional coupling between mechanical fields and material structure, in which the initial mechanical field drove density changes, density changes modified the local elastic modulus and changes in modulus in turn altered the mechanical field. This feedback loop progressively reorganised the internal distribution of stress and strain. Comparison between the initial and final configurations showed preservation of the overall spatial organisation of mechanical domains, redistribution of stress intensity in response to local stiffness changes and attenuation of extreme local variations through iterative adaptation. The final configuration maintained a gradient-based structure, although with modified intensity and distribution. No uniform or global increase or decrease in density was observed; rather, the model consistently produced localised, gradient-driven changes characterised by heterogeneous adaptation patterns.

Notably, regions sharing similar absolute σh values but differing in local |∇σ_h| magnitude exhibited markedly different density evolution: high-gradient zones underwent measurable Δρ while mechanically homogeneous zones with comparable σ_h remained essentially unchanged. This dissociation between absolute magnitude and adaptive response supports the interpretation that, within the present formulation, structural evolution depends on the spatial context of the mechanical field rather than on its global magnitude. Taken together, the iterative simulations demonstrate that bone density evolution in this model is governed by spatial gradients of mechanical variables; that the system evolves towards a mechanically balanced, non-uniform configuration; that adaptation occurs in a localised manner rather than as a homogeneous global response; and that global force magnitude alone does not determine the resulting remodelling pattern. The spatial distribution of hydrostatic stress, the corresponding gradient vector field and the harmonic deformation of the gingival envelope are illustrated in Fig. 3.

Figure 3: Computational validation of the Python post-processing pipeline. (A) Hydrostatic stress σ_h (RdBu_r: blue compression, red tension) ranging from −0.0767 MPa (cervical compression) to +0.1007 MPa (opposing wall tension); green line marks σ_h = 0. (B) Gradient vector field ∇σ_h (arrows) with log-scaled magnitude background (hot); mean |∇σ_h| = 0.0061 MPa/mm. (C) Harmonic mesh deformation of the gingival envelope after 30 increments, showing smooth adaptation without distortion. Scale: mm.

These results extend the static gradient-based framework described in Part I into a time-dependent mechanobiological context, showing how spatially organised mechanical fields may influence the evolution of bone properties under sustained orthodontic loading [8].

Pipeline Validation and Computational Benchmark

The Python post-processing pipeline successfully ingested the discrete geometry and stress fields transferred from the ANSYS Workbench environment. The extraction, spatial indexing and gradient computation modules executed with 100% numerical stability across all benchmark cycles, with no convergence failures or mesh corruption events.

In a control sample of 500 cortical elements oriented toward vestibular tipping, the software identified a physical hydrostatic stress range from −0.0767 MPa (maximum focal compression at the cervical/crestal margin) to +0.1007 MPa (maximum tension on the opposing alveolar wall). These values are consistent with the low-to-moderate stress regime characteristic of aligner-mediated adult orthodontic loading reported in Part I [8]. The KD-Tree robustly resolved element adjacency within the 1.5 mm search radius, with an average query time of less than 0.1 s per node on a standard workstation configuration.

The spatial derivative computation yielded a mean mechanical gradient magnitude of 0.0061 MPa/mm. This local gradient vector provided the three-dimensional rate-of-change and trajectory required by the geometric engine. The harmonic deformation linear system assimilated these gradient-driven displacement vectors without registering mesh distortions, folds or inverted elements, enabling continuous, smooth and anatomically coherent morphing of the gingival envelope adapted to the new dental positioning. Fig. 3 illustrates the spatial distribution of hydrostatic stress and the corresponding gradient vector field over the cortical domain for a representative tipping increment. The convergence behaviour of the iterative remodelling algorithm, confirming numerical stability and mechanical redistribution without oscillatory instability, is shown in Fig. 4.

Figure 4: Convergence behaviour of the iterative remodelling algorithm. (A) Maximum density update |Δρ|_max per iteration, showing monotonic decay toward a quasi-stationary plateau (0.002 g/cm³). (B) Peak |σ_h|_max (blue) and |∇σ_h|_max (orange) across 30 iterations, both attenuating as local stiffness redistribution modifies load transfer. Absence of oscillatory instability confirms numerical robustness.

Taken together, these benchmark data confirm that the closed-loop FEA-Python pipeline is computationally robust and geometrically stable. The spatial gradient magnitudes obtained (∼0.006 MPa/mm) are of the same order of magnitude as the pressure differentials reported in poromicromechanical analyses of osteocyte lacunar flow, supporting the biological plausibility of the gradient-based stimulus descriptor Φ employed in the iterative model [16]. The successful harmonic morphing of 30 consecutive 0.25-mm aligner increments demonstrates that the pipeline satisfies the prerequisites for integration into patient-specific orthodontic planning software (Part III).

Discussion

The present study extends the framework developed in Part I from a spatial description of the mechanical field to its temporal evolution under sustained orthodontic loading [8]. While the static analysis identified strain and stress gradients as defining features of the mechanical environment in adult alveolar bone, the iterative model implemented here demonstrates how these gradients may drive the progressive redistribution of tissue properties over successive simulation cycles. The results indicate that, within the adopted formulation, structural adaptation is not organised according to absolute values of stress or strain, but according to the spatial structure of the mechanical field, which remains the dominant feature throughout the iterative process. This observation is consistent with the gradient-based interpretation of orthodontic biomechanics proposed in Part I and with previous experimental and computational evidence indicating that bone adaptation in adult contexts is governed by the heterogeneity and directional structure of the mechanical environment rather than by force magnitude alone [2,5,6,8,20,29].

Classical mechanostat-based interpretations of bone adaptation assume well-defined strain thresholds and a relatively direct, deterministic coupling between stimulus magnitude and biological response [1]. However, the behaviour observed in the present iterative model does not conform to such assumptions. No abrupt threshold-driven transitions were observed; adaptation occurred in a gradual and localised manner; and regions exhibiting similar stress or strain magnitudes produced different responses depending on their spatial context. This behaviour is consistent with experimental observations in adult orthodontic bone, where adaptive capacity is reduced, biological response is variable and no direct proportionality exists between sustained force magnitude and remodelling outcome [5,6,23-25]. The age-related attenuation of bone metabolism in adult patients, the reduced osteocyte mechanosensitivity associated with deterioration of the lacuno-canalicular network and the diminished responsiveness of periodontal ligament cells have all been described as factors contributing to a constrained adaptive capacity in mature alveolar bone [23-27]. Within this context, mechanostat-based formulations developed primarily on the basis of long-bone remodelling may not adequately describe the behaviour of adult alveolar cortical bone under sustained orthodontic loading.

A consistent finding of the present model was the spatial correlation between hydrostatic stress gradients (|∇σ_h|) and local density changes (Δρ), suggesting that the relevant mechanical stimulus in this formulation is not only the local stress state but its spatial variation. From a mechanobiological perspective, this observation is compatible with current evidence regarding osteocyte mechanosensing. Macroscopic mechanical loading generates lacunar pressure gradients within the lacuno-canalicular network, which in turn drive interstitial fluid flow at nanoscale velocities around osteocytes [16,17]. Osteocytes are highly sensitive to fluid shear stress and to local pressure differentials and respond through the release of paracrine factors that modulate osteoblast and osteoclast activity [5,13,14,15]. The spatial organisation of the macroscopic mechanical field therefore acts as an organising principle of the local mechanobiological environment: regions of elevated gradient magnitude correspond to regions where lacunar-canalicular flow is most likely to be perturbed, defining zones of differential adaptive potential. This interpretation is consistent with poromicromechanical analyses showing that physiological bone strains generate undrained lacunar pressures sufficient to stimulate osteocyte mechanosensing [16]. Importantly, the presence of mechanically plausible gradients does not imply guaranteed remodelling, particularly in adult alveolar bone, where biological response remains conditional and patient-specific [2,3,8].

At the molecular level, recent advances in mechanobiology provide additional support for the relevance of gradient-based mechanical descriptors. Mechanosensitive ion channels such as PIEZO1 have been identified in periodontal ligament cells, osteocytes and bone cells and have been implicated in the conversion of mechanical stimuli into biochemical signals during orthodontic tooth movement [3,4,18,19]. Activation of PIEZO1 in periodontal ligament cells has been shown to promote alveolar bone remodelling on the tension side of orthodontic loading and to modulate osteogenic and osteoclastic activity [4,18]. The role of PIEZO1 in osteocytes during orthodontic tooth movement is currently under investigation, with conflicting evidence regarding the extent to which osteocyte-specific deletion alters tooth movement in animal models [7]. While the present finite element formulation does not explicitly model these molecular pathways, the spatial organisation of hydrostatic stress and its gradients provides a mechanically plausible substrate for such mechanotransduction processes to operate. The consistency between the spatial structure of the mechanical field and the regions of differential density evolution supports a framework in which mechanical gradients delineate domains of potential biological sensitivity, within which the activation of cellular mechanotransduction pathways is conditional on local biological context.

The iterative framework employed in this study exhibits a key behavioural feature: bidirectional coupling between mechanics and structure. The initial mechanical field drives density changes, which in turn modify the local elastic modulus, which alters the subsequent mechanical field. This feedback loop progressively reorganises the internal distribution of stress and strain. Rather than converging towards a uniform state, the system evolves towards a mechanically structured configuration in which gradients persist while extreme local variations are attenuated. This behaviour is broadly consistent with iterative orthodontic finite element models published in the literature, which have similarly reported gradual evolution of tooth-supporting tissues under sustained loading [9,10,11,21]. However, in those formulations the driving variable is typically strain energy density or PDL strain and the relevance of spatial gradients within cortical bone is not explicitly addressed. The framework adopted here differs in that the mechanical stimulus function Φ explicitly incorporates the spatial gradient of hydrostatic stress as a primary descriptor, in line with continuum remodelling models that distinguish hydrostatic and deviatoric components of mechanical state [7,13,20,29]. This formulation does not aim to replace strain-based models but to complement them by providing an alternative framework better suited to the constrained, near-homeostatic context of adult alveolar cortical bone, where strain magnitudes typically remain within low-to-moderate ranges and the most informative aspect of the mechanical environment is its spatial organisation rather than its absolute intensity [2,6,7,8].

From a biomechanical and clinical perspective, these results suggest that the nature of orthodontic loading and its spatial distribution may be more relevant than force magnitude alone in shaping adult alveolar bone response. Loading configurations that minimise extreme gradients while maintaining controlled tooth movement may favour mechanical stability over large redistributions of tissue properties. In this sense, strategies such as root-controlled tipping or moment-dominant loading [8], modulation of force systems through aligner staging [11] or careful selection of attachments and force vectors [10,28,29] can be interpreted as methods to shape the spatial gradient field rather than simply increasing or decreasing applied forces. This interpretation aligns with the conditional and spatially constrained view of Wolff’s law proposed in Part I and reinforces the recommendation to consider the structural anatomy of adult cortical bone – particularly the thin buccal and lingual cortical plates – when designing orthodontic mechanics in adult patients [8,26,29].

A methodological feature of the present work, inherited from Part I, is the explicit separation between mechanical computation and biological interpretation [8]. Finite element analysis was used solely to compute mechanical fields, while the definition of remodelling-relevant variables and the iterative density update were performed through external post-processing. This separation avoids embedding speculative biological laws within the numerical solver and is consistent with recommendations to keep mechanical and biological assumptions distinct when modelling complex adaptive systems [2,13]. The model does not propose a new biological law, but provides a controlled mechanical framework for analysing the temporal interaction between spatially organised mechanical stimuli and simplified structural adaptation in adult alveolar cortical bone.

The present work is subject to several limitations. The finite element simulations rely on homogeneous, isotropic and linearly elastic assumptions for the bone compartments, which do not capture the anisotropic osteon architecture and regional mineralisation gradients of adult cortical bone. Although the periodontal ligament was modelled with a generalised viscoelastic formulation (α_i^G = 0.925, α_i^K = 0.890, τ_i = 0.15 s), the present implementation does not incorporate full non-linear hyperelasticity, fibre-reinforced anisotropy or spatially variable PDL thickness, which may influence long-term stress relaxation, load transfer asymmetry and the temporal evolution of hydrostatic stress gradients under sustained orthodontic activation. The remodelling parameters (k, n, Δt) were selected to ensure numerical stability rather than calibrated against experimental or clinical data. No explicit biological time scale, cellular mechanism or biochemical process is included in the model and no histological, micro-CT or longitudinal clinical validation has been performed. Accordingly, the results should be interpreted as a mechanobiological exploration of how spatially organised mechanical stimuli may interact with simplified structural evolution, rather than as a predictive simulation of clinical outcomes. Future work should address the sensitivity of the model to parameter variation, integrate experimental data on adult alveolar bone density distributions [23-26] and incorporate more refined descriptions of the periodontal ligament and the lacuno-canalicular network when computational resources allow [16,17]. Furthermore, a native limitation of the current iterative framework is the assumption of a constant alignment displacement vector across all 30 cycles. In clinical scenarios, an aesthetic aligner progressively loses its elastic deflection capacity as the target tooth displaces. Consequently, the 30 simulated iterations presented herein should be interpreted as the local, short-term mechanobiological evolution of the mature alveolar tissue within a single activation stage (0.25 mm), isolating the initial tissue remodeling dynamics prior to macro-scale tooth kinematics.

This study represents the second step of a structured research programme. Part I established the spatial organisation of the mechanical field under orthodontic loading; the present work extends that description into a time-dependent mechanobiological model based on iterative finite element simulation; and a subsequent study (Part III) will translate this framework into orthodontic movement planning software, allowing patient-specific simulation of treatment progression with explicit consideration of mechanical gradients [8]. Within this sequence, the present model functions as a bridge between mechanical description and computational application. The results support a framework in which adult alveolar bone adaptation under orthodontic loading is governed by mechanical gradients rather than by absolute stimulus magnitudes; structural evolution emerges from iterative mechanical-material coupling rather than from threshold-driven responses; and the system tends to converge towards organised, non-uniform configurations rather than homogeneous redistributions of tissue properties. From this perspective, Wolff’s law in adult orthodontic conditions can be reinterpreted not as a deterministic relationship between force magnitude and bone remodelling, but as a gradient-driven mechanical framework that defines regions of potential adaptation within a constrained and patient-specific biological context [1,2,8].

To translate the computational stability verified in this iterative study into clinician-focused diagnostic sequences, a standalone geometric pre-processing and visualization framework was developed. As a preliminary architectural validation, a dedicated spatial rendering module was implemented to compute and display real-time root-to-bone proximity envelopes and initial boundary gradient approximations across both maxillary and mandibular arches (Fig. 5). This clean front-end visualization structure handles the spatial K-Dimensional Tree (KD-Tree) neighbor indexing required to analyze complex patient-specific anatomical boundary conditions prior to deploying the definitive, non-linear physical calculations within an external structural finite element solver (CalculiX), effectively establishing the operational groundwork for full patient-specific multi-tooth orthodontic staging in Part III.

Figure 5: Translational framework for multi-arch spatial visualization and boundary pre-processing (Prerequisite for Part III). Three-dimensional rendering of patient-specific maxillary and mandibular structures demonstrating the clinical orientation interface. The color-coded cortical maps represent preliminary geometric distance and proximity gradient approximations between the dental roots and the surrounding alveolar cortical bone envelopes. This visualization core operates as a spatial pre-processing engine, structuring the complex multi-tooth boundary conditions and local proximity thresholds before transmitting the anatomical dataset to an external non-linear finite element core (CalculiX) for definitive mechanobiological simulation.

Conclusion

An iterative finite element framework integrating mechanical computation with external numerical updating was successfully implemented to simulate the temporal evolution of bone density and elastic modulus in adult alveolar cortical bone under sustained orthodontic loading. The proposed model extends the static, gradient-based mechanical framework established in Part I into a time-dependent mechanobiological formulation. Within this formulation, structural evolution is governed primarily by the spatial organisation of the mechanical field – particularly by hydrostatic stress gradients rather than by absolute magnitudes of stress or strain. The system evolves towards stable, non-uniform configurations characterised by organised mechanical domains and progressive attenuation of extreme local variations, without exhibiting threshold-driven or deterministic remodelling behaviour. These findings reinforce a non-deterministic, context-dependent interpretation of bone adaptation in adult orthodontic conditions and support a gradient-based reinterpretation of Wolff’s law in this clinical setting. The present work does not aim to predict clinical outcomes; rather, it provides a mechanically consistent framework for the subsequent integration of mechanobiological reasoning into orthodontic treatment planning software, to be developed in Part III.

Conflict of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding Statement

This research did not receive any specific grant from funding agencies in the public, commercial or non-profit sectors.

Acknowledgement

The authors have no acknowledgments to declare.

Data Availability Statement

The data supporting the findings of this study are available from the corresponding author upon reasonable request.

Ethical Statement

The project did not meet the definition of human subject research under the preview of the IRB according to federal regulations and therefore was exempt.

Informed Consent Statement

Not applicable.

Authors’ Contributions

All authors contributed equally to this paper.

References

1. Frost HM. A 2003 update of bone physiology and Wolff’s law for clinicians. Angle Orthod. 2004;74(1):3-15.
2. Li Y, Jacox LA, Little SH, Ko CC. Biomechanical and biological responses of periodontium in orthodontic tooth movement: Update in a new decade. Int J Oral Sci. 2021;13(1):20.
3. Jiang Y, Liu X, Chen L, Wang Y, Li Z, Zhang H, et al. Mechanosensitive Piezo1 in periodontal ligament cells promotes alveolar bone remodeling during orthodontic tooth movement. Front Physiol. 2021;12:767136.
4. Jiang Y, Lin H, Chen Y, Lan Y, Wang H, Li T, et al. Piezo1 contributes to alveolar bone remodeling by activating β-catenin under compressive stress. Am J Orthod Dentofacial Orthop. 2024;165(4):458-70.
5. Klein-Nulend J, Bakker AD, Bacabac RG, Vatsa A, Weinbaum S. Mechanosensation and transduction in osteocytes. Bone. 2013;54(2):182-90.
6. Wang J, Huang Y, Chen F, Li W. The age-related effects on orthodontic tooth movement and the surrounding periodontal environment. Front Physiol. 2024;15:1460168.
7. Nottmeier C, Lavicky J, Lopez MG, Knauth S, Kahl-Nieke B, Amling M, et al. Mechanical-induced bone remodeling does not depend on Piezo1 in dentoalveolar hard tissue. Sci Rep. 2023;13:9563.
8. Salas G. Mechanical gradients in adult alveolar cortical bone under orthodontic tipping: a finite element study (I). J Dent Health Oral Res. 2026;7(1):1-13.
9. Wang Y, Gao L, Liu Y. Numerical simulation of mandible bone remodeling under tooth loading: A parametric study. Sci Rep. 2019;9:14887.
10. Wu JL, Liu YF, Zhuang JR, Liu CK, Jiang CP. Evaluation of long-term orthodontic tooth movement considering bone remodeling process and in the presence of alveolar bone loss using finite element method. Orthod Waves. 2016;75(4):75-84.
11. Liu Y, Hu W. The effect of maxillary molar distalization with clear aligner: A 4D finite-element study with staging simulation. Prog Orthod. 2023;24(1):14.
12. Horina J, Wagner D, Cherian J, Cevidanes L, Steiniger B, Wriggers P, et al. An in-silico approach to modeling orthodontic tooth movement using stimulus-induced external bone adaptation. J Mech Behav Biomed Mater. 2021;124:104891.
13. Branecka N, Yildizdag ME, Ciallella A, Giorgio I. Bone remodeling process based on hydrostatic and deviatoric strain mechano-sensing. Biomimetics. 2022;7(2):59.
14. Qin L, Liu W, Cao H, Xiao G. Molecular mechanosensors in osteocytes. Bone Res. 2020;8:23.
15. Bonewald LF. The amazing osteocyte. J Bone Miner Res. 2011;26(2):229-38.
16. Scheiner S, Pivonka P, Hellmich C. Poromicromechanics reveals that physiological bone strains induce osteocyte-stimulating lacunar pressure. Biomech Model Mechanobiol. 2016;15(1):9-28.
17. Lin C, Karim L, Wang J, et al. Mechanical deformations of bone generate interstitial fluid flow at nanoscale velocities around osteocytes. Front Bioeng Biotechnol. 2025;13:1639788.
18. Du Y, Yang K. Role of mechanosensitive ion channel Piezo1 in tension-side orthodontic alveolar bone remodeling in rats. Arch Oral Biol. 2023;155:105798.
19. Zhou T, Gao B, Fan Y, Liu Y, Feng S, Cong Q, et al. Piezo1/2 mediate mechanotransduction essential for bone formation through concerted activation of NFAT-YAP1-β-catenin. eLife. 2020;9:e52779.
20. Gross TS, Edwards JL, McLeod KJ, Rubin CT. Strain gradients correlate with sites of periosteal bone formation. J Bone Miner Res. 1997;12(10):1603-13.
21. Hamanaka R, Yamaoka S, Anh TN, Tominaga JY, Koga Y, Yoshida N. Numerical simulation of orthodontic bone remodeling. Orthod Waves. 2017;76(4):189-96.
22. Schulte FA, Ruffoni D, Lambers FM, Christen D, Webster DJ, Kuhn G, et al. Strain-adaptive in silico modeling of bone adaptation: A computer simulation validated by in-vivo micro-computed tomography data. PLoS One. 2013;8(11):e79522.
23. Sun Y, Huang J, Yang Z, Su X, Wang B, Hu Y, et al. Age-related alveolar bone maladaptation in adult orthodontics: finding new ways out. Int J Oral Sci. 2024;16(1):25.
24. Xu J, Huang Z, Liu W, Jiang J, Liang J, Wu C, et al. Changes in alveolar bone structure during orthodontic tooth movement in adolescent and adult rats: A microcomputed tomography study. Orthod Craniofac Res. 2023;26(3):494-503.
25. Zheng Y, Zhu C, Zhu M, Lei L. Difference in the alveolar bone remodeling between the adolescents and adults during upper incisor retraction: A retrospective study. Sci Rep. 2022;12:9161.
26. Ren Y, Maltha JC, Kuijpers-Jagtman AM. Age-dependent responses of alveolar bone to orthodontic force in rats. Am J Orthod Dentofacial Orthop. 2019;156(5):620-30.
27. Danz JC, Degen M. Selective modulation of the bone remodeling regulatory system through orthodontic tooth movement — a review. Front Oral Health. 2025;6:1472711.
28. Cattaneo PM, Dalstra M, Melsen B. The finite element method: A tool to study orthodontic tooth movement. J Dent Res. 2005;84(5):428-33.
29. Tiwari AK, Goyal A, Prasad J. Modeling cortical bone adaptation using strain gradients. Proc Inst Mech Eng H. 2021;235(5):544-54.