Abstract
Fluoride-containing formulations are widely used in preventive dentistry; however, whether clinically used formulations exhibit distinct protein-level interaction tendencies that could be relevant to enamel biology remains insufficiently explored. In this study, we comparatively assessed four fluoride-containing formulations (sodium fluoride, amine fluoride/olaflur, sodium monofluorophosphate, and sodium difluorophosphate) against enamel-associated enzymes and matrix proteins using an in silico workflow. Molecular docking (AutoDock Vina) was performed against carbonic anhydrase II, tissue-nonspecific alkaline phosphatase, and matrix metalloproteinase-20, along with AlphaFold models of key enamel matrix proteins for hypothesis generation. Across targets, olaflur showed the most favorable docking scores and recurrent stabilizing contacts. Based on the docking rank and biological relevance, the olaflur– matrix metalloproteinase-20 complex was subjected to a 100-ns molecular dynamics simulation, demonstrating stable backbone deviation, limited residue-level fluctuations near the binding region, persistent hydrogen bonding, and favorable end-state binding free energy by Molecular Mechanics/Generalized Born Surface Area (MM/GBSA). Collectively, these exploratory in silico results generate hypotheses suggesting that organic amine fluorides may display distinct molecular-level interaction tendencies with enamel-associated proteins, providing a basis for future experimental studies to examine their potential biological relevance.
Introduction
Dental caries is one of the most common chronic diseases in children worldwide. While there are differences by region, socioeconomic factors, diet, oral hygiene, and preventive care access, dental caries continues to be extremely common and burdensome. It warrants investigation of preventive measures. 1 The widespread distribution of caries, along with its early appearance and multifactorial causes, creates long-term health effects, which collectively demonstrate its significant burden on public health. 2 Children are particularly susceptible to caries because of established cariogenic biofilms, high exposure to fermentable carbohydrates, and immature enamel in early life. 3 Consequently, preventive strategies for enamel protection and early lesion arrest are the mainstays of pediatric dental care. The use of fluoride-based therapies is the most widely recommended preventive intervention for controlling early enamel demineralization in children. 4 Fluoride contributes to enamel remineralization, reduces mineral loss, and interferes with bacterial metabolism in dental biofilms. Fluoride serves as a fundamental element in pediatric dental treatment when used topically in various professional and non-professional forms, such as toothpaste, varnishes, gels, mouth rinses, and fluoride-releasing dental materials and sealants.5,6 Available fluoride-based agents include sodium fluoride, sodium monofluorophosphate, sodium difluorophosphate, and amine fluorides (e.g. olaflur). Although they share a similar preventive goal, these compounds vary in terms of chemical structure, dissociation and oral bioavailability, enamel substantivity, and antibacterial efficacy, which may contribute to differences in their clinical performance. Most fluoride research uses mineral-focused mechanisms to interpret results without adequately examining the interactions affected by specific formulations. 7
Research has primarily investigated how fluoride interacts physicochemically with enamel mineral phases through fluorapatite development, calcium–phosphate precipitation control, and decreased enamel solubility at acidic levels. 8 Although these mineral-centric mechanisms have been well established, enamel is not a purely inorganic material. The highly ordered crystal morphology results from a highly regulated biomineralization process that occurs during amelogenesis and is directed by enamel matrix proteins (EMPs). 9 EMPs, including amelogenin, ameloblastin, and enamelin, control crystal nucleation, spatial orientation, and structural maturation. In a complementary process, matrix metalloproteinase-20 (MMP20) serves as the initial protease involved in the processing of amelogenin during the secretory stage of enamel formation and has been identified as critical for normal enamel formation and matrix remodeling. 9 Some enzymes that are not classified as matrix proteins play crucial roles in maintaining enamel homeostasis. Carbonic anhydrase II (CA2) regulates local pH via bicarbonate buffering. The activity directly impacts enamel demineralization and remineralization, while CA2 was included as a biologically relevant protein to study the enamel microenvironment, considering the pH/ionic dependency of fluoride behavior in solution. Tissue-nonspecific alkaline phosphatase (ALPL) plays a key role in regulating mineral deposition by controlling local phosphate concentrations via pyrophosphate hydrolysis. 10 The role of this enzyme in enamel and bone mineral maturation has been well described and supports the inclusion of F-releasing agents that show differences in their interaction behavior with enzymatic regulators of mineralization.
During enamel maturation, most EMPs undergo degradation and removal; however, small protein fragments are retained within the interprismatic enamel and persist throughout the lifetime of the tooth. 11 These residual components are more prominent in primary teeth and newly erupted permanent enamel, giving pediatric enamel biological characteristics that differ from those of fully matured adult enamel. Proteomic and structural studies have indicated that retained enamel proteins can influence enamel permeability, remineralization behavior, and responses to externally applied chemical agents.12,13 Limited research exists on the direct interaction of fluoride-releasing agents with EMPs. Experimental studies have shown that fluoride exposure influences protein expression profiles, protein degradation rates, and mineral–protein interactions during enamel formation. These effects may contribute to differences in remineralization outcomes, variability in protein processing following fluoride treatment, and, in certain contexts, to developmental disturbances such as enamel hypomineralization or fluorosis. Differences among fluoride formulations, including amine and phosphate-based fluorides, may reflect agent-specific interactions with enamel proteins that extend beyond conventional mineral-based mechanisms.
Recent developments in structural bioinformatics and computational chemistry have enabled the examination of interactions between small molecules and proteins at the molecular scale. Approaches such as molecular docking allow the comparison of binding affinity, preferred interaction regions, and stabilization of non-covalent contacts, while related computational methods can be used to examine interaction stability under defined conditions. 14 Similar multiscale computational strategies have increasingly been applied in chemically complex biomaterial and biointerface research to generate mechanistic hypotheses prior to experimental validation.15,16 Simultaneously, the increasing availability of experimentally resolved protein structures and high-confidence predicted models, together with curated chemical structures from public databases, has enabled systematic in silico analysis of biologically relevant protein–ligand systems. These approaches are particularly applicable to pediatric enamel, where residual matrix proteins are more abundant and fluoride exposure occurs repeatedly during preventive care procedures. In addition to the major EMPs, several enzymes involved in enamel homeostasis and mineral regulation are also present and functionally relevant. Matrix metalloproteinase-20 (MMP20), CA2, and tissue-nonspecific ALPL contribute to protein processing, local pH regulation, and phosphate availability in the enamel environment.17 –19 Although these proteins are central to enamel development and maturation, their interactions with clinically used fluoride-releasing agents have not been comparatively examined.
Despite these advancements and the clear relevance of protein interactions, a comparative evaluation of formulation-dependent fluoride interaction tendencies with enamel-associated proteins remains limited, particularly for enzymatic regulators that shape enamel mineral chemistry and matrix remodeling. Therefore, to address this gap, we applied a standardized in silico docking workflow across key enamel-associated proteins, followed by molecular dynamics simulation and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) analysis of the top-ranked complex to generate mechanistic hypotheses on how different fluoride formulations may differ at the protein interface in enamel biology.
From a chemical research perspective, this study explores the effects of formulation-based variances in fluoride speciation, structure, polarity, and functional groups on non-covalent binding affinities to biologically relevant protein surfaces, placing the research at the nexus of computational chemistry, biomolecular interactions, and preventive dentistry.
Materials and methods
This study followed a stepwise in silico workflow comprising comparative molecular docking of clinically used fluoride formulations against enamel-associated proteins, followed by molecular dynamics simulation and MM/GBSA analysis of the top-ranked complex to evaluate the interaction stability and energetic favorability.
Selection and retrieval of protein structures
The three-dimensional structures of enamel-associated proteins were selected based on their biological relevance to enamel formation, maturation, and mineral regulation. The experimentally resolved protein structures were retrieved from the RCSB Protein Data Bank (https://www.rcsb.org/). MMP20 (PDB ID: 2JSD), CA2 (PDB ID: 3KS3), and tissue-nonspecific ALPL (PDB ID: 7YIW) were selected due to their established roles in enamel matrix processing, pH regulation, and phosphate metabolism, respectively. 20 When experimentally determined structures were unavailable, high-confidence predicted models for amelogenin (AMELX), enamelin (ENAM), and ameloblastin (AMBN) were obtained from the AlphaFold database (https://alphafold.ebi.ac.uk/). AlphaFold models serve only to generate hypotheses and facilitate exploratory research. Since the models lack experimentally derived ligand-binding sites and may include errors in surface topology and local side-chain orientation, caution was used when interpreting the results from AlphaFold models and they were not considered to be equally confident as results from experimentally determined structures. 21
Protein preparation
The Discovery Studio Visualizer (version 2024) was used for protein preparation. All protein structures were checked for completeness and integrity before docking. Removal of water molecules, co-crystallized ligands, non-potential heteroatoms, and ions that were not involved in maintaining the structural integrity was performed to ensure that there was no possibility of interfering with the ligand-binding site. The prepared protein structures were saved in PDB file format and converted into PDBQT format using PyRx (version 0.8) software. 22 During this process, polar hydrogen atoms were added and Gasteiger partial charges were assigned. All protein targets were prepared in the same manner for consistency during the docking process.
Ligand selection and preparation
Four fluoride-releasing materials commonly used in pediatric and preventive dentistry were selected as ligands. The ligand 3D structures were downloaded directly from PubChem in SDF format. The ligands were prepared using PyRx (version 0.8) for energy minimization purposes. The Universal Force Field (UFF) was used as the force field using a conjugate gradient method to achieve an root-mean-square deviation (RMSD) of <0.001 kcal (mol·Å)−1 within 500 steps. 22 After energy minimization, the ligand file format was converted to PDBQT, and the rotatable bonds were assigned. The docking process allows researchers to evaluate ligands against each other under identical experimental conditions. Docking of NaF was performed solely as a computational approximation for the presence of fluoride-induced electrostatic contact during scoring rather than as an estimation of a putative docked complex. As NaF dissociates into solvated ions, no direct comparison between the docking scores of NaF and undissociated ligands (organic or phosphate) should be made. The comparison should be qualitative.
Binding site identification
The research identified experimentally bound pockets through the analysis of protein structures obtained from experimental crystal determination. The PrankWeb server predicted putative binding pockets from experimentally determined crystal structures by combining structural data with machine learning–based features to identify potential binding cavities. 22 The highest probability pocket for each experimentally determined structure was chosen, and its center coordinates were used to center the docking grid. The absence of experimentally confirmed binding sites for AlphaFold-predicted protein structures required the docking grids to extend across the entire predicted protein structure to ensure an unbiased ligand interaction analysis across its surface area. The docking grid box sizes were determined by extending the grid dimensions to include the entire predicted binding pocket. The grid sizes were between 22 and 28 Å along each axis, depending on the predicted pocket shape. The docking grid box sizes for the AlphaFold-predicted structures were selected to cover the entire volume of the protein structure. Grid boxes were centered on the centroid of each protein and extended by 5 Å in each direction beyond the maximum extent of the protein structure, allowing for an unbiased sampling of potential interaction regions throughout the full protein surface area. For AlphaFold-predicted EMPs, whole-protein docking was applied in an exploratory manner and should not be considered a replacement for pocket-guided docking. The choice to perform whole-protein docking instead of attempting to pinpoint potential binding sites arose from the need to avoid overinterpreting predictions without experimental evidence to support the presence of these sites. The docking performance of AlphaFold structural models functions as a qualitative measure of surface interaction capability but remains incomparable to binding scores derived from other structural systems.
Molecular docking simulation and interaction analysis
Docking simulations were performed using AutoDock Vina on the PyRx platform.23,24 A docking exhaustiveness of 20 was used to strike a reasonable balance between thorough conformational sampling and computational expense for all protein–ligand complexes. Because of the approximate scoring in AutoDock Vina that can differ with search exhaustiveness, docking was performed in this work for the purpose of relative ranking and comparison of interaction patterns rather than for discriminating closely spaced scores. Small score differences especially those within the margin of error were therefore treated with skepticism and not considered significant evidence of true energetic discrimination without further validation. For each protein–ligand combination, multiple binding poses were generated, and the pose associated with the most favorable predicted interaction energy was retained for comparative evaluation. The docked complexes were subsequently examined using the Discovery Studio Visualizer. Hydrogen bonding, electrostatic interactions, and other noncovalent contacts were identified based on established geometric and distance-based criteria. Interaction patterns were compared across different protein targets and fluoride-releasing agents to evaluate the relative binding behavior of enamel-associated proteins. 24
Molecular dynamics simulation
Molecular dynamics simulations were performed to investigate the structural stability and dynamic behavior of the Olaflur-2JSD (MMP20) complex using GROMACS 2025.1. This complex was selected based on its comparatively favorable docking profile and biological relevance to the evaluated protein–ligand systems. Accordingly, the molecular dynamics simulation was intentionally restricted to a single representative complex, and no inference regarding the dynamic stability of other fluoride–protein interactions was intended within the scope of this study. As the analysis was limited to a single 100-ns run without replicates, the MD data should be understood as providing supportive evidence for the behavior of this specific complex, rather than definitive confirmation of complete dynamic convergence. The protein was modeled using the CHARMM36 force field, and the Olaflur parameters were generated using SwissParam. The protonation states corresponding to pH 7.4 were assigned using PDB2PQR, after which the complex was solvated in a CHARMM-TIP3P water box with a minimum 10 Å buffer. The system was neutralized and adjusted to an ionic strength of 0.15 M using Na+/Cl− ions to mimic physiological conditions. 24 Energy minimization was performed using the steepest descent algorithm until the maximum force converged to below 1000 kJ mol⁻¹ nm⁻¹. Subsequent equilibration and production runs were conducted according to standard molecular dynamics protocols.
Equilibration was conducted under the constant number of particles, volume, and temperatur (NVT) ensemble for 500 ps at 300 K using the v-rescale thermostat, 25 followed by constant number of particles, pressure, and temperature (NPT) equilibration for 500 ps at 1 bar using the Parrinello–Rahman barostat. 26 Periodic boundary conditions were applied throughout, with long-range electrostatics treated using the Particle Mesh Ewald method and a 1.2 nm cut-off for short-range interactions. Production MD was performed for 100 ns using a 2 fs time step, and the trajectories were saved every 100 ps. Structural stability and flexibility were evaluated using RMSD, root-mean-square fluctuation (RMSF), radius of gyration, solvent-accessible surface area, and hydrogen-bond analyses with standard GROMACS tools. 27
All structural analyses (RMSD, RMSF, Rg, solvent-accessible surface area (SASA), and hydrogen bonds) were computed using the same production trajectory window specified in the Results section to maintain consistency. Unless otherwise stated, the analyses summarize the full 100-ns production run.27,28
MM/GBSA binding free energy analysis
Binding free energy estimates of the Olaflur–2JSD complex were calculated using the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) method as implemented in gmx_MMPBSA v1.6.3, which combines GROMACS trajectories with AmberTools for end-state free energy estimation.29,30 The total binding free energy (ΔG_bind) was decomposed into molecular mechanical energy components, including van der Waals (ΔE_vdW) and electrostatic (ΔE_ele) interactions, along with solvation free energy terms comprising polar (ΔG_polar) and nonpolar (ΔG_nonpolar) contributions, according to the following equation:
The polar solvation energy was computed using the Generalized Born implicit solvent model with GB-Neck2 parameterization (igb = 5), whereas the nonpolar solvation term was estimated from the solvent-accessible surface area using the Linear Combination of Pairwise Overlaps (LCPO) method. MM/GBSA calculations were not incorporated with entropic terms since the current study was performed to provide a relative end-state prediction of binding energetics rather than being used to calculate an absolute free-energy value. This limitation is particularly relevant for olaflur because its flexible structure may incur a substantial configurational entropy penalty upon binding, which could make the apparent favorability of the calculated ΔG_bind more negative than the true thermodynamic free energy. The computed ΔG_bind values are predominantly enthalpic in nature and are only to be considered as qualitative representations of the binding affinity. Considerable cancelation of the positive electrostatic and negative polar solvation terms is expected for ligand–protein complexes involving polar interactions. MM/GBSA calculations in these systems should not be viewed at the component level, and the total value will depend on the implicit solvent model/dielectric assumptions (Genheden & Ryde, 2015). A total of 1001 snapshots were extracted at uniform intervals from the 100 ns-equilibrated production trajectory of each protein–ligand complex. MM/GBSA calculations were performed separately for the complex, receptor, and ligand, and the final binding free energy was obtained as:
Results
Enamel-associated protein structures and fluoride ligands
An overview of the enamel-associated protein targets and fluoride-containing ligands included in this study is provided in Tables 1 and 2, respectively. The protein set included experimentally determined structures of MMP20, CA2, and tissue-nonspecific ALPL, along with AlphaFold-predicted models of amelogenin, enamelin, and ameloblastin. Considering both enzymatic proteins and structural matrix components allowed for a comparison of fluoride interaction behavior across proteins involved in different stages of enamel formation and maturation.
Structural information on the enamel-associated proteins used in this study.
Note. Experimentally resolved structures were retrieved from the RCSB PDB, and AlphaFold-predicted models were obtained from the AlphaFold database (version 6 for human proteins).
Fluoride-containing compounds selected for docking analysis.
The ligand set included sodium fluoride, Olaflur, sodium difluorophosphate, and sodium monofluorophosphate, which represent clinically relevant fluoride formulations with distinct chemical properties. The selection of these compounds allowed for a systematic assessment of how differences in molecular structure and functional groups influence their interaction behavior with enamel-associated proteins.
Predicted binding sites of enzymatic protein targets
The predicted ligand-binding pockets identified in the experimentally resolved protein structures are summarized in Table 3. In MMP20 (PDB ID: 2JSD), the highest-ranked pocket was associated with a pocket score of 27.80 and a probability score of 0.903. This pocket comprised a continuous stretch of residues spanning A_180 to A_249 and was localized to a compact region with center coordinates at (4.0930, −4.5180, 6.6950). For CA2 (PDB ID: 3KS3), the top-ranked binding site exhibited the highest pocket score among the evaluated enzymes (31.24), with a probability of 0.923. The predicted pocket included residues distributed between positions A_92 and A_209, and the calculated pocket center was located at (−4.8518, 2.9112, 14.3102). For tissue-nonspecific ALPL (PDB ID: 7YIW), the predicted binding pocket displayed a pocket score of 17.62 and a probability of 0.803. Residues contributing to this pocket were distributed across chains C and D, with the pocket center positioned at (23.2897, 103.8318, 87.9746). These sites were used as reference regions in subsequent docking analyses. For EMPs represented by AlphaFold-predicted models, including amelogenin, enamelin, and ameloblastin, predefined ligand-binding pockets were unavailable owing to the absence of experimentally validated binding information. Docking calculations were performed using grids that encompassed the full protein structures, allowing ligand interactions to be explored across the entire protein surface.
Predicted active site pockets using PrankWeb server.
Docking binding energies of fluoride-releasing agents
The docking binding energies for the complete set of protein–ligand complexes are presented in Table 4. Significant variations in the predicted interaction tendencies were observed among the evaluated fluoride-containing compounds. For the enzymatic protein targets, Olaflur consistently produced the most favorable predicted energies, with values of −5.4 kcal mol−1 for tissue-nonspecific ALPL, −5.2 kcal mol−1 for CA2, and −5.1 kcal mol−1 for MMP20, as shown in Table 4. These values represent the most favorable docking scores obtained across all ligand–protein combinations examined. These values represent the most favorable docking scores obtained across all ligand–protein combinations examined. However, the range of the enzymatic Olaflur docking scores was narrow (−5.1 to −5.4 kcal mol−1) and should not be overinterpreted as indicating a true energetic hierarchy, given the approximate nature and inherent uncertainty of the docking scoring functions. These results were interpreted as indicating a broadly favorable interaction tendency across the evaluated enzymatic targets rather than a strong quantitative ranking among them. Phosphate-containing fluorides exhibited intermediate binding. The sodium monofluorophosphate and sodium difluorophosphate energies ranged from approximately −3.3 to −4.7 kcal mol−1 (exact values varied by protein), whereas the energies for sodium fluoride were consistently high (i.e. unfavorable) for each protein at approximately −2.0 kcal mol−1 or above. These values are not representative of actual molecular binding (as discussed in the Methods section) but rather a reflection of a computational proxy for the transient electrostatic affinity of the fluoride ion and serve as a methodological baseline. Sodium fluoride belongs to a distinct chemical class from the other ligands assessed in this study. Sodium fluoride readily dissociates into Na+ and F− ions in aqueous solution and does not exist as a molecule with stable directional bonding. Consequently, the docking results obtained for sodium fluoride should not be interpreted as evidence of stable protein–ligand complex formation in vivo. Instead, they reflect a computational approximation of the transient electrostatic contact propensity of the fluoride ion under the applied scoring function. These values are not directly comparable to the docking scores obtained for structurally intact organic or phosphate-based fluoride compounds and are included here only as a qualitative reference to highlight the methodological limitations of applying conventional docking frameworks to freely solvated ions. As such, NaF numbers should not be compared outside the scope of the same protocol.
Predicted docking binding energies of fluoride-containing agents with enamel-associated proteins. The docking values for sodium fluoride reflect a simplified computational approximation of the fluoride ion interaction propensity and do not represent stable protein–ligand binding. These values are not directly comparable to the docking scores obtained for intact molecular ligands, such as Olaflur or phosphate-based fluorides.
Note. The results obtained for AlphaFold-predicted enamel matrix proteins are exploratory and reflect generalized surface interaction tendencies rather than site-specific binding predictions. These values are presented separately for hypothesis generation and should not be interpreted with the same confidence as the results derived from experimentally resolved protein structures.
For EMPs represented by AlphaFold-predicted models, broadly similar directional patterns were observed, with Olaflur generally yielding more favorable docking scores than inorganic and phosphate-based fluorides. However, these observations should be interpreted as qualitative and exploratory, as the absence of experimentally validated binding pockets and the use of whole-protein docking limit the specificity and confidence of the interaction predictions derived from these models. Thus, data collected from AlphaFold models of EMPs should be treated as low-confidence hypothesis-generating data, rather than direct evidence equivalent to experimentally determined protein structures.
Residue-level interaction analysis of Olaflur
The residue-level interaction profiles of olaflur with selected enzymatic protein targets are summarized in Table 5. A representative three-dimensional docking conformation of the olaflur–MMP20 (2JSD) complex is shown in Figure 1. Across all complexes, the interactions involved a combination of hydrogen bonding and additional non-covalent contacts with polar and hydrophobic residues located within the predicted binding regions.
Key intermolecular hydrogen-bonding and non-covalent interactions of Olaflur with selected protein targets.
Note. Hydrogen bonding and non-covalent interactions were assigned during post-docking analysis in Discovery Studio Visualizer using standard geometric and distance-based criteria. Non-covalent interactions include hydrophobic contacts (alkyl and π–alkyl interactions), π–π stacking, π–cation interactions, electrostatic contacts, and van der Waals interactions.
Representative docking conformations and interaction profiles of Olaflur with selected enzymatic protein targets: (a) matrix metalloproteinase-20 (MMP20; PDB ID: 2JSD), (b) carbonic anhydrase II (CA2; PDB ID: 3KS3), and (c) tissue-nonspecific alkaline phosphatase (ALPL; PDB ID: 7YIW).
In the MMP20–olaflur complex (PDB ID: 2JSD), hydrogen bonds were observed with GLU277, THR247, TYR248, and VAL223. These contacts were accompanied by non-covalent interactions involving TYR180, ALA194, HIS236, PHE193, and HIS191. The selected residues established a local interaction space with polar and aromatic amino acids, which supported Olaflur in maintaining its orientation inside the predicted pocket. For CA2 (PDB ID: 3KS3), hydrogen bonds were formed between GLN92 and HIS119. Additional non-covalent contacts were observed with ILE91 and PHE131, residues that neighbor the hydrogen-bonding residues. These interactions positioned Olaflur in a localized area consisting of hydrophobic and polar residues. Hydrogen bonds were found between ASP337 and ASP294 in the tissue-nonspecific ALPL (PDB ID: 7YIW) complex. The stabilization of the docked ligand was supported through non-covalent interactions with MET401, LEU449, HIS451, HIS341, and TYR338. The spread of interacting residues throughout this region shows how Olaflur settles on the protein surface instead of forming a single-point contact.
RMSD
At the beginning of the simulation, the apo protein (represented by the black curve) displayed an initial RMSD increase, which was followed by significant changes during the trajectory, maintaining RMSD values mostly between ~0.32 and 0.45 nm, as shown in Figure 2. This indicates that the protein has higher conformational flexibility in the absence of the ligand. The Olaflur–2JSD complex (red curve) rapidly equilibrated in the first few nanoseconds of the simulation, and the RMSD values remained stable with small fluctuations. The overall average RMSD of the complex was 0.38 nm, with a maximum deviation of ~0.48 nm. This implies a fairly stable trajectory behavior for the representative docked complex in a given time window. The smaller fluctuation amplitude in the complex compared to the apo protein suggests that Olaflur binding may contribute to the structural stabilization of 2JSD under simulated conditions.
Backbone RMSD profiles of apo 2JSD (black) and the Olaflur-2JSD complex (red).
RMSF
The RMSF plot analysis of residue mobility within the Olaflur-2JSD complex is shown in Figure 3. RMSF analysis of the Olaflur–2JSD complex revealed generally low fluctuation at the residue level, with the vast majority of values ranging between 0.05 and 0.15 nm and an average RMSF of 0.13 nm. Minor localized peaks with higher values (~0.55 nm) are also noted, which suggests that certain small regions of the protein did not lose flexibility during the simulation. Since residue-by-residue annotation of these maxima was not performed as part of the present analysis, these peaks can only be conservatively interpreted as being attributed to flexible peripheral residues/protein unfolding rather than global protein instability. Notably, no systematic increase in RMSF was detected throughout the predicted site of ligand interaction.
RMSF plot of the Olaflur-2JSD complex.
Rg and SASA analysis
The Rg values were relatively constant along the trajectory, with an average value of 1.59 nm, suggesting the preservation of the overall folded state (Figure 4(a)). The observed plot variations throughout the simulation, measured in nanometers, represent typical protein breathing behavior rather than large-scale unfolding. The flat Rg profile suggests that olaflur binding does not perturb the global architecture of the 2JSD protein and is consistent with the formation of a compact and stable complex. The SASA profile showed moderate fluctuations along the trajectory, with values ranging from 88.99 to 105.04 nm² and an average value of 96.93 nm² (Figure 4(b)). The stability of SASA throughout the simulation demonstrated that ligand attachment did not cause significant structural changes or protein unfolding.
Structural compactness and solvent exposure of the Olaflur-2JSD complex during molecular dynamics simulation: (a) radius of gyration (Rg) profile and (b) solvent-accessible surface area (SASA) profile.
Hydrogen bond analysis
Hydrogen bonding between Olaflur and 2JSD complex was examined over time to see how long these interactions persisted during the 100 ns trajectory (Figure 5). The number of hydrogen bonds between the ligand and the protein fluctuated mostly between 2 and 5 over the entire trajectory, and these fluctuations were short in nature and could be due to the thermal motion of the entire complex. Over the course of the trajectory, the ligand was never observed to completely dissociate from the binding region. This suggests that the representative complex remained associated throughout the simulation. Hydrogen bonds that form the repeating pattern shown here help maintain the docked pose throughout the simulation time frame but are not conclusive evidence of bound stability outside the current simulation.
Hydrogen bond dynamics of the Olaflur-2JSD complex.
MM/GBSA binding free energy analysis
The MM/GBSA binding free energy components of the olaflur–2JSD complex throughout the equilibrated MD trajectory were analyzed to determine the energetic factors influencing ligand binding (Figure 5). The predicted binding free energy was −18.14 kcal mol−1, indicating that the end-state interaction between Olaflur and 2JSD was enthalpically favorable in the MM/GBSA model. The favorable net value resulted primarily from significant cancelation between a large favorable electrostatic term and a large unfavorable polar solvation term, which is typical of MM/GBSA assessments of polar ligand–protein complexes. Because of this phenomenon, care should be taken when over-interpreting individual energy terms in isolation. Furthermore, the endpoint ΔG_bind estimate is affected by the implicit-solvent model and dielectric assumptions employed. For these reasons, MM/GBSA results are best considered as qualitative, enthalpy-dominated endpoint predictions of interaction favorability, not quantitative thermodynamic binding free energies. This value also does not include entropic contributions, and for a flexible ligand such as olaflur, the loss of configurational entropy upon binding may be substantial; therefore, the reported favorable ΔG_bind is likely overestimated relative to the true thermodynamic free energy. It should be emphasized that the MM/GBSA results reported herein are specific to the Olaflur–MMP20 complex and pertain to the energetic behavior of this system under dynamic conditions. Further MD simulations are required before making assumptions about the stability or energetic favorability of other complexes or fluoride formulations.
Binding energy decomposition indicated that the electrostatic (ΔE_ele = −268.23 kcal mol−1) and van der Waals (ΔE_vdW = −20.17 kcal mol−1) interactions are the dominant intermolecular forces responsible for binding, which produced an overall, strongly favorable gas-phase interaction energy (ΔG_gas = −288.40 kcal mol−1). The desolvation of the charged species upon binding is responsible for the polar solvation term (ΔG_GB = +274.08 kcal mol−1), while the nonpolar term (ΔG_surf = −3.82 kcal mol−1) was a small, stabilizing factor in the overall solvation free energy (ΔG_solv = +270.26 kcal mol−1). A negative ΔGbind value indicates favorable end-state interaction energetics for the olaflur–2JSD complex within the MM/GBSA calculations.
Discussion
This study aimed to conduct a molecular-level comparison between clinically significant fluoride-releasing substances and essential enamel-related proteins involved in enamel development and mineral regulation. While traditional measures of fluoride efficacy have been related to mineral-centric processes, this study computationally examined whether alternative fluoride species may demonstrate unique binding affinities at the protein level in enamel systems. The biological makeup of pediatric enamel is of particular interest because it maintains a higher quantity of enamel residual matrix proteins.
Of the compounds studied, Olaflur had the most favorable predicted docking scores within the confines of the applied silico model; however, this should not be viewed directly as significant proof of increased biological or clinical efficacy. As demonstrated in the docking results, olaflur had the most favorable predicted docking scores with experimentally solved enzyme targets. AlphaFold-predicted models for EMPs demonstrated tentative similar interaction propensities, but these results are speculative without verification and depend heavily on the current models and whole-protein docking methodology. This observation is consistent with the chemical properties of amine fluorides, which have an organic backbone that can establish multiple hydrogen bonding interactions and stabilize non-covalent interactions with the protein surface.
The interaction patterns discovered for EMPs with AlphaFold-based models need to be regarded with less confidence than patterns recognized from experimentally established enzyme targets. The lack of identified binding sites required the docking of the full proteins, which is more likely to result in non-specific surface contacts. Models originating from AlphaFold data should be used as preliminary indications to help form hypotheses regarding preferential binding mechanisms.
Unlike organic and phosphate-based fluorides, sodium fluoride (NaF) fails to maintain its molecular structure when it enters a biological environment. In water or other aqueous solutions, NaF quickly dissociates to yield fluoride anions (F−), the active species that helps remineralize enamel and combat bacteria. Traditional docking algorithms, which evaluate stable, conformationally rigid ligands, do not adequately function within the applicability domain that includes freely solvated ions because of their molecular docking limitations. The docking score for NaF is expected to mostly reflect weak, non-directional electrostatic interactions. The NaF results demonstrate that ions cannot be directly compared using docking, while establishing them as methodological controls instead of suggesting inferior binding to molecular fluoride forms.
Sodium monofluorophosphate and sodium difluorophosphate displayed intermediate interaction behaviors between the fluoride ion and Olaflur. The former have a larger molecular structure and contain functional groups for hydrogen and electrostatic bonding interactions, which might be the reason for their relatively high docking scores. The data suggest that the molecular size and availability of such functional groups influence protein association behavior under docking conditions.
A possible reason for this difference is that olaflur has an amphiphilic organic scaffold that can present several weak interactions simultaneously (hydrogen bonding, electrostatics, and hydrophobic contacts), which may lead to a longer residence time in the shallow surface pockets than in the small inorganic/ionic species. Simple fluoride salts experience solvation and ionic dissociation in water-based solutions, leading to poor representation in classical docking scoring, which accounts for their weak and non-directional interaction patterns found in in silico studies.
Residue-level interactions provided additional support for these observations. Olaflur contacted polar, charged, and hydrophobic residues in the putative binding sites of MMP20, CA2, and ALPL, respectively. In the MMP20 complex, contacts with GLU277, THR247, and TYR248 located the compound in a region linked to protein-processing activity. For CA2, interactions with GLN92 and HIS119 positioned the ligand near residues involved in enzymatic and buffering activity, whereas the contacts observed in the ALPL complex involved residues related to phosphate coordination and structural stabilization. These interaction patterns suggest that Olaflur can adapt to diverse protein environments related to enamel biology under the docking conditions applied.
To better understand the dynamic features of the Olaflur–MMP20 complex in physiological environments, we conducted molecular dynamics simulations of the complex. The lower RMSD fluctuations of the complex compared to the apo protein can be interpreted as a stabilizing effect of the ligand. The RMSF plot (Figure 6) shows that the flexibility at the residue level is minimal. The predicted binding region exhibited low fluctuations during the simulation, whereas the loop and terminal regions exhibited high fluctuations. The binding mode remained constant throughout the simulation, demonstrating that the interaction was not transient.
MM/GBSA binding energy components of the olaflur–2JSD complex.
MM/GBSA calculation for the olaflur–MMP20 complex yielded a negative value for the free binding energy of the complex, which indicates a favorable energetic profile under the chosen computational conditions. The calculated interaction was governed by electrostatic and van der Waals forces that overcame the unfavorable entropic polar solvation energy loss (i.e. desolvation penalty) of the charged groups on the receptor upon complex formation. MM/GBSA calculations do not determine exact free binding energies but the results suggest energetic viability for the stable formation of Olaflur–MMP20 complexes within the simulation.
The molecular dynamics and MM/GBSA simulations focused only on the Olaflur–MMP20 complex to carry out a thorough representative analysis rather than a broad comparative dynamics research. Although this complex stability and positive energetics indicate the likelihood of this interaction being maintained, these results do not imply the potential of other Olaflur–protein pairs and other fluoride formulations tested in this study. Only one top complex was used for the molecular dynamics simulations and MM/GBSA calculations, and only one 100-ns production run was utilized for the dynamic analysis. Therefore, dynamic observations should not be extrapolated to all fluoride–protein complexes. Conclusions regarding relative dynamic stabilities would necessitate multiple simulations and analyses of other ligand–protein complexes.
The present data reveal that numerous fluoride species demonstrate specific interaction capabilities with enamel-related proteins at the molecular scale, based on our current computational model. This suggests potential ways in which proteins might bind, although it does not confirm the direct biological effects. Experimental verification of the predicted modes of interaction is required to confirm whether this binding affects enzyme activity, protein stability, or other local biochemical conditions pertinent to enamel. Binding assays, spectroscopic studies, biochemical validation, and structural and clinically relevant studies are needed to confidently ascertain biological significance. The increased presence of residual matrix proteins in children’s enamel suggests that they could play a key role in the response of enamel to repeated fluoride exposure during dental prevention programs.31,32
The results of this study should be interpreted in the context of its limitations. Since this was a purely computational study, the findings require further confirmation via laboratory work, animal testing, and clinical trials (if applicable) before we can draw any biological or preventive conclusions. AlphaFold-predicted structures were used for the EMPs because of the absence of experimentally resolved structures of these proteins at full length, which may have affected the binding-site predictions. Docking scores represent relative affinities for interaction but not absolute binding free energies. Docking scores also do not consider the effects of physiological variables such as pH, salivary components, biofilms, and the constantly changing enamel surface. Physiological conditions cause rapid dissociation of inorganic fluoride salts, and the subsequent use of fluoride ions as stand-ins during docking studies prevents accurate modeling of ionic solvation impacts.
AlphaFold predictions for EMPs require cautious interpretation because of the inherent model uncertainty combined with the use of unbiased docking methods. These should not be equated with findings originating from experimentally determined structures.
MM/GBSA analysis does not consider configurational and conformational entropy, which influences binding free energy calculations mainly for flexible ligands and surface-associated binding modes. This issue is especially important in the present study because olaflur is a flexible ligand, and omission of entropy may materially exaggerate the apparent energetic favorability of the olaflur–MMP20 interaction. Therefore, the ΔG_bind values are best thought of as a measure of enthalpy-dominated tendencies toward interaction and are likely to overestimate the absolute binding strength. Normal mode analysis paired with quasi-harmonic entropy estimation or alchemical free energy approaches will enable future research to achieve a more detailed thermodynamic analysis of fluoride–protein interactions.
Running molecular dynamics simulations exclusively on one top-ranked complex reduces the dynamic comparisons between various ligands and protein targets and restricts the stability and MM/GBSA conclusions to the Olaflur–MMP20 system.
This study presents a structure-based approach to evaluate protein interactions with various fluoride formulations under specified constraints. By including both modulatory enzymes and EMPs, our findings suggest that protein–fluoride interactions should be considered along with mineral-phase effects to understand fluoride behavior in enamel. In addition, a direct comparison between molecular ligands and free solvated species, such as fluoride ions, is limited within the classical docking approach and should be considered only on a qualitative basis. This study establishes a foundation for future experimental work to confirm protein interactions that depend on the formulation and to understand their significance in pediatric enamel biology.
From a preventive dentistry perspective, formulation-dependent protein interactions observed in silico could provide a molecular basis for further investigation into the reported differences among fluoride-containing formulations rather than serving as direct explanations for clinical outcomes. Such molecular-level associations may, in theory, influence protein-mediated processes at the enamel surface, including enzymatic regulation or matrix interactions, which can be evaluated experimentally in future studies. However, the present study did not directly model pediatric enamel physiology, and any relevance to pediatric contexts remains speculative and requires biological validation.
Conclusion
Through an exploratory computational workflow, we assessed the propensity of a panel of clinically available fluoride formulations to interact with enamel-related proteins. Olaflur achieved the most optimal docking scores against the assessed targets; however, the scores of the experimentally resolved enzymatic targets were similar (−5.1 to −5.4 kcal mol−1) and are too closely ranged to infer strong energetic trends. Overall, the chosen Olaflur–MMP20 complex was stable throughout the 100 ns simulation (average backbone RMSD = 0.38 nm) and displayed a favorable end-point MM/GBSA estimate (ΔG_bind = −18.14 kcal mol−1). As exploratory endpoints, these results do not demonstrate biochemical, biological, or clinical precedence of any one fluoride formulation. They do however suggest hypotheses to be explored experimentally to assess whether formulation-dependent fluoride–protein interactions translate to measurably relevant effects in enamel systems.
Footnotes
Acknowledgements
All the authors are thankful to the King Khalid University, Saudi Arabia, for the financial Support. The authors would like to acknowledge Mr S. Sadhu for his valuable technical assistance in conducting the molecular docking simulations and for providing support in ligand and protein preparation workflows.
Ethical considerations
Not applicable.
Consent to participate
Not applicable.
Consent for publication
Not applicable.
Author contributions
Funding
The authors extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University for funding this work through Small Research Project under grant number RGP1/161 /46.
Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Data availability statement
The data supporting the findings of this study are available from the corresponding author upon reasonable request.