Multiscale Modeling of Bovine Chymosin Catalysis: From AlphaFold3 and Replicated Molecular Dynamics to QM/MM Characterization of the Pre-Reactive State

1. Introduction

Chymosin (EC 3.4.23.4) is an aspartic protease expressed in the stomach of young mammals (class Mammalia) and belongs to the family of pepsin-like aspartic proteases [1,2,3].

The primary biological substrate for chymosin (Chn) is one of the main milk proteins: kappa-casein (κ-CN). In this molecule, the enzyme attacks selectively and rapidly a single Chn-sensitive or key bond (in bovine κ-CN, this is the Phe105-Met106 bond). Hydrolysis of the Chn-sensitive bond initiates the milk coagulation process and leads to the formation of a milk curd. The transformation of milk into a curd slows the movement of food from the stomach to the intestine, promotes gentle stimulation of the newborn digestive tract, and ensures the complete absorption of the complex of milk nutrients [1,4,5,6].

Due to their narrow substrate specificity, chymosins stand out from other gastric pepsin-like proteases because of their high milk-clotting activity (MCA) and relatively low overall proteolytic activity (PA). The dominance of MCA over total PA determines the special biological role of these enzymes: Chns not only participate in digestive processes in newborn mammals but also contribute to the development of their immune defenses [4].

Most placental animals are characterized by postnatal transfer of maternal antibodies (Ab) to the newborn via milk [7,8,9]. Synthesis of Chn by cells of the gastric mucosa begins in the late prenatal period. A special biological role of Chns is to create a weak proteolytic background in the neonatal stomach, which protects immunoglobulins and other protective proteins of milk from degradation. As a result, intact antibodies bypass the stomach and enter the intestine, where they are transcellularly transferred into the blood and lymph directly, providing immune protection to the newborn [7,8].

Preparations of proteolytic enzymes have been used by humans for thousands of years to initiate the first stage of cheese making: milk coagulation and obtaining a milk curd suitable for further technological processing. Historically, most milk-clotting enzyme preparations were extracted from the stomachs of ruminant farm animals (cows, sheep, goats) and received the general name: rennets (R). The main active R preparations are chymosin and pepsin (EC 3.4.23.1). The most widely used preparations are bovine Rs [1,9,10].

In addition to a high MCA/PA ratio, which ensures preferential hydrolysis of the key peptide bond in κ-CN, chymosins are characterized by technologically attractive indices of thermal stability (TS) and dependence of coagulation activity on Ca²⁺ concentration and pH. This combination of unique biochemical properties has allowed Chns to surpass many other milk-coagulating proteases and become the primary milk-clotting enzyme for cheese making [1,9,11].

Chymosin is synthesized as an inactive precursor, preprochymosin, consisting of 381 amino acid (a.a.) residues with a molecular weight (Mm) of ~43.0 kDa. The “pre” fragment, or signal peptide, consists of 16 aa and exhibits hydrophobic properties. During translation, the signal peptide passes through the membrane of the endoplasmic reticulum (ER), “pulling” a molecule of preprochymosin. In the ER, the signal peptide is eliminated and preprochymosin is converted into an inactive zymogen, prochymosin (365 aa, Mm of ~40.8 kDa). The “pro” fragment, or activation peptide (42 aa), blocks the active site of the enzyme, which creates a steric hindrance to substrate binding. From the ER, prochymosin is transported to the Golgi apparatus and is subsequently secreted from the cell. The formation of active Chn (323 aa, Mm of ~36.0 kDa) occurs at acidic pH values due to autocatalytic hydrolysis, during which the activation peptide is removed from the N-terminus of prochymosin [12,13,14].

There are two known genetic versions of Chn: A and B, which differ by one amino acid residue at position 244: Asp244 (Chn A) → Gly244 (Chn B) [15]. Chn B predominates in R preparations. The Asp244Gly substitution leads to a change in the physicochemical properties of the enzyme: in contrast to Chn A, variant B is less active, has increased TS and is more stable at pH < 3.5 [16]. The pH optima of chymosins A and B differ and are, respectively, 4.2 and 3.7 [17]. The term “chymosin” used in this paper, unless otherwise stated, refers to bovine Chn B (Bos taurus).

Several crystal structures of bovine chymosin have been determined by X-ray diffraction [18,19,20,21]. In addition, the Protein Data Bank (PDB) includes the three-dimensional structure of Chn in complex with the inhibitor norstatin, also determined by X-ray diffraction (PDB ID: 1CZI) [22].

No experimentally determined three-dimensional structure of full-length bovine κ-CN is currently available in the PDB, which may be related to the rheomorphic nature of this protein [5,15,23].

Chymosin is a monomeric globular protein consisting of 323 amino acid residues [14]. The residue numbering used in this study follows that of the X-ray structure deposited as PDB ID: 1CMS [24]. The spatial organization of Chn is characteristic of aspartic proteases, among which a high degree of structural homology is observed [25,26]. Chymosin has a pseudosymmetric bi-lobal structure with a fold known as a “β-barrel” [27]. The enzyme molecule (dimensions: ≈40 × 60 × 65 Å) consists of two domains, the N- and C-terminal domains, separated by a cavity (or cleft) approximately 2.5 Å wide [18]. The cavity contains a substrate-binding site and two catalytic aspartate residues, Asp34 and Asp216, which are oriented toward each other (“face-to-face”) and form part of the conserved Asp-Thr-Gly motifs (Figure 1). The position of catalytic aspartates is stabilized by a network of hydrogen bonds, where the key role is played by two Thr residues, whose interaction is called “the fireman’s grip” [19,25]. The enzyme contains three disulfide bridges (Cys47–Cys52, Cys207–Cys211, and Cys250–Cys283) and a cis-proline residue (Pro25) [16].

In the immediate vicinity of the active site of Chn, there is a region of the polypeptide chain 71-82, forming a β-hairpin, which carries a Tyr residue at position 77. According to X-ray diffraction data, depending on the conformation of the R-group of Tyr77, the enzyme can exist in two alternative forms: open (active) (PDB ID: 3CMS) or closed (self-inhibited) (PDB ID: 1CMS, 3CMS, 4CMS) [29]. In the self-inhibited form (Chn-SI), the side group of Tyr77 occupies part of the binding pocket, where the substrate is located during interaction with the enzyme. In the active form (Chn-Act), the side chain of Tyr77 is “folded back” and the substrate-binding site of the enzyme is free [29,30]. It is assumed that the P8-P4 fragment of κ-CN with the His98-Pro99-His100-Pro101-His102 structure, called the histidine-proline (His-Pro) cluster, can act as an allosteric effector that stimulates the conversion of the closed form of Chn into the open form [17]. Earlier kinetic analysis of bovine κ-CN-derived tryptic fragments spanning residues 98–111 and 98–112 showed that inclusion of the His-Pro-rich 98–102 segment markedly improved overall substrate properties toward chymosin [31]. According to [30], preliminary incubation of Chn with a pentapeptide repeating the structure of the His-Pro cluster increases the proteolytic activity of the enzyme in relation to the κ-CN (f103-109) peptide by approximately 189 times.

A number of works deal with the study of the mechanism of catalytic cleavage of the amide bond by aspartic proteases [2,3,18,19,20,21,24,25,26,27,32,33].

It has been proposed that hydrolysis of the scissile peptide bond proceeds through nucleophilic attack on the peptide carbonyl carbon by a water molecule [2,34,35,36]. However, there is no consensus on the precise details of the hydrolytic mechanism [2,17,29,32].

For example, the authors of [37] proposed the reaction mechanism in which the carbonyl bond (C=O) of Asp34, together with the hydrolyzable bond Phe105–Met106, are protonated and simultaneously subjected to nucleophilic attack by a water molecule, which is activated by deprotonated Asp216 (Figure 2).

According to this proposed mechanism, the hydrolysis process begins with proton transfer from a water molecule to the Asp34/Asp216 catalytic dyad (Figure 2, structure I). Simultaneously, a proton migrates from the dyad to carbonyl oxygen of Phe105 of the cleaved peptide bond Phe105–Met106 (Figure 2, structure II). Proton migration is the driving force for the bond formation between the substrate and water and leads to the formation of intermediate tetrahedral oxyanion (II) [37,38]. The resulting oxyanion is hydrolyzed by an acid-base catalysis mechanism, with simultaneous proton transfer from the intermediate compound (II) to the aspartate dyad and from this dyad to the nitrogen atom of the substrate, forming hydrolysis products (III). A similar mechanism has been previously described in [24,28,29,30].

The role of water molecules in the structure and functioning of aspartic proteases is also discussed in [39,40].

The model of κ-CN and Chn interaction is studied in [41], using molecular docking, conformational search algorithms and molecular dynamics (MD). This model suggests that Lys111 and Lys112 residues of κ-CN interact with the N-terminal domain of Chn, replacing a water molecule, and the histidine-proline cluster of the substrate binds to the C-terminal domain of the enzyme.

If the assumption that the allosteric activator of Chn is a fragment of homologous κ-CN is correct, then it provides an idea of why the bovine enzyme does not always efficiently hydrolyze heterologous κ-CNs containing amino acid residues other than the His-Pro cluster [42,43]. Although the process of allosteric activation of Chn is widely discussed in the literature [2,19,28,29,30,41,44], this mechanism has not been elucidated at the molecular level due to the complexity of its detailed description.

The role of the N- and C-terminal sequences of the Chn-sensitive region of κ-CN (residues 96–119) [45], including the His-Pro cluster, in the proposed regulation of chymosin activity remains incompletely understood. In particular, it is unclear how the open and self-inhibited conformations of Chn are populated under solution conditions and to what extent substrate binding contributes to stabilization of the open state. The molecular organization of the Chn/κ-CN complex immediately preceding cleavage of the Phe105–Met106 bond also remains insufficiently characterized. In particular, the positioning and mobility of the scissile region relative to the Asp34/Asp216 catalytic dyad, as well as the local organization of water molecules within the catalytic microenvironment, require further atomistic characterization. A detailed description of these structural features may help to clarify the molecular basis of chymosin substrate recognition and specificity.

In the present study, we applied a multiscale computational strategy combining AlphaFold3-based structural modeling, molecular docking, replicated molecular dynamics simulations, and QM/MM calculations to characterize the pre-reactive Chn/κ-CN complex. Particular attention was paid to the conformational behavior of Tyr77, the contribution of the His-Pro cluster to stabilization of the enzyme–substrate complex, the positioning and mobility of the Phe105–Met106 region relative to the catalytic dyad, and the organization of the hydrogen-bonded water network within the active site. QM/MM refinement was used to examine the local catalytic microenvironment, while an exploratory one-dimensional relaxed scan along the O(water)···C(carbonyl) coordinate was employed to probe the energetic response of the active-site environment to progressive water approach toward the scissile carbonyl group. The resulting multiscale model provides a structural framework for interpreting the high specificity of bovine chymosin toward its homologous substrate and may also contribute to understanding differences in its interactions with heterologous κ-caseins [42,43].

2. Materials and Methods

2.1. Protein and Peptide Preparation

Three starting models of bovine chymosin were considered: the X-ray structures representing the self-inhibited and open conformations (PDB IDs: 4AA8 and 4AUC, respectively) [46,47], and a structure predicted using AlphaFold3 (Chn-Fold) [48]. The X-ray structures differ in the conformation of the flexible loop comprising residues 72–79, which contains Tyr77 and is associated with the self-inhibited (Chn-SI) and open (Chn-Act) states of the enzyme.

The experimental chymosin structures were prepared for subsequent calculations by adding the missing Asn291–His292–Ser293 residues, reverting the Val111Phe substitution to the native residue [19,49], defining the Cys47–Cys52, Cys207–Cys211, and Cys250–Cys283 disulfide bonds, adding hydrogen atoms, removing crystallographic water molecules, and performing restrained structural optimization using the OPLS3e force field at pH 6.5 [50]. The Chn-Fold model was subjected to the same protonation and structural-preparation protocol where applicable.

The amino acid sequence of bovine κ-CN (NP776719.1) was obtained from the GenBank database. The Chn-sensitive region consisting of 16 aa was chosen for the study (residues 97–112) [51]. To predict the three-dimensional structure of the selected fragment, the AlphaFold3 methodology was used, implemented in the web application (https://alphafoldserver.com) [48]. The predicted local distance difference test (pLDDT) score provided by AlphaFold3 was used as a per-residue confidence metric for the predicted structure. Higher pLDDT values indicate greater confidence in the local geometry of the predicted model [48,52].

2.2. Molecular Docking

Molecular docking of the Chn-sensitive κ-CN fragment was performed using the representative structure of the most populated cluster obtained from the selected Chn-Fold MD-10 trajectory as the receptor. The κ-CN peptide model generated with AlphaFold3 was used as the ligand. Docking calculations were carried out using the PIPER protein–protein docking protocol [53]. During pose generation, the receptor and peptide were treated as rigid bodies, and the ligand was systematically rotated and translated relative to the receptor. The resulting docking solutions were evaluated using the PIPER scoring function and their geometric compatibility with the chymosin catalytic cleft.

A positional constraint was applied to retain the Phe105–Met106 region of κ-CN within 7 Å of the Asp34/Asp216 catalytic dyad [1,51]. Docking poses satisfying this geometric criterion were subsequently examined for spatial accommodation of the peptide within the binding cleft and for the relative orientation of the scissile Phe105–Met106 region with respect to the catalytic residues. The selected docking pose was used as the starting structure for subsequent molecular dynamics simulations of the Chn-Fold + κ-CN complex.

2.3. Molecular Dynamics Simulations

Molecular dynamics (MD) simulations were performed for three models of bovine chymosin: Chn-Act, Chn-SI, and Chn-Fold. In addition, the selected Chn-Fold + κ-CN docking complex was subjected to MD simulations to characterize the stability and conformational dynamics of the enzyme–substrate complex. All simulations were performed using the Desmond molecular dynamics engine within the Maestro environment [54,55].

Each system was placed in an explicit-solvent box with a minimum buffer distance of 15 Å between the solute and the box boundaries. The systems were solvated using the TIP3P water model and supplemented with NaCl to a concentration of 0.15 M. Simulations were carried out in the NPT ensemble at 310 K and 1.01325 bar using the OPLS3e force field [50]. Protein structures were prepared at pH 6.5 prior to MD simulations.

For each modeled system, 10 independent 300-ns MD trajectories were generated. Each replicate was initialized independently using randomized initial atomic velocities generated by the simulation software. Thus, the total sampling time for each system was 3.0 μs.

For the Chn-Fold system, the time-averaged RMSD and the corresponding standard deviation were calculated separately for each of the 10 independent trajectories. The MD-10 trajectory, which showed the lowest mean RMSD (2.123 Å) together with moderate conformational fluctuations (σ = 0.19 Å), was selected for subsequent structural clustering. RMSD-based clustering was performed using the 300-ns MD-10 trajectory only. For each cluster, the representative frame was defined as the structure with the lowest average RMSD relative to the other members of that cluster. The representative structure of the most populated MD-10 cluster was subsequently used as the receptor for molecular docking of the κ-CN peptide.

For the Chn-Fold + κ-CN complex, RMSD characteristics were likewise evaluated separately for the 10 independent trajectories. The MD-5 trajectory was selected for subsequent structural analysis based on its RMSD behavior and conformational stability. Analysis of enzyme–substrate interactions and RMSD-based structural clustering were performed using the selected 300-ns MD-5 trajectory only. The representative structure of the most populated MD-5 cluster was used for structural characterization of the pre-reactive complex.

2.4. Quantum Mechanics/Molecular Mechanics (QM/MM) Calculations

A frame from the selected 300-ns MD-5 trajectory of the Chn-Fold + κ-CN complex was chosen as the starting structure for QM/MM calculations. The frame was selected based on the open-state geometry of Tyr77 and a short distance between the Phe105 scissile region and Asp34. The local catalytic region was subsequently refined using a hybrid QM/MM approach.

The QM region comprised the catalytic residues Asp34 and Asp216, the Phe105–Met106 region of κ-CN, and four water molecules located in the catalytic microenvironment. The QM subsystem was described using the B3LYP hybrid density functional [56,57,58] with the 6-31G(d,p) basis set, whereas the remaining part of the system was treated at the MM level using the OPLS2005 force field [59].

To characterize the distribution of water molecules near the catalytic region, water occupancy was analyzed in the volume surrounding Phe105, Met106, Asp34, and Asp216 using the Volume Map plugin of VMD [60]. A water molecule located within a region of high solvent occupancy was selected as the putative nucleophilic water for the QM/MM model.

An exploratory one-dimensional relaxed QM/MM scan was then performed along the distance between the oxygen atom of the selected water molecule and the carbonyl carbon of Phe105, O(water)···C(carbonyl). At each point of the scan, the remaining degrees of freedom were allowed to relax according to the applied QM/MM optimization protocol. The resulting energy profile was used to evaluate the energetic response of the modeled active-site environment to progressive water approach. The scan was not intended to determine a complete reaction pathway or to locate a transition state.

3. Results

3.1. Analysis of Chymosin and Peptide Structures

It is assumed that Chn can exist in two conformations (forms): self-inhibited or active ones [28,29]. In the self-inhibited form, the Tyr77 residue covers the active site of the enzyme (torsion angle θ_{N-Cα-Cβ-Cγ} = 175°) and hinders access to it by the Chn-sensitive sequence of κ-CN, while in the active form, on the contrary, it does not hinder the interaction of catalytic aspartates with the substrate (angle θ = 300°) (Figure 3).

Mature κ-CN contains 169 aa [62]. The Chn-sensitive region of bovine κ-CN, consisting of 16 a.a., which carries a histidine-proline cluster and the key residues Phe105 and Met106, was chosen as a substrate model [63] (Figure 4).

AlphaFold3 was used to predict the three-dimensional structure of the Chn-sensitive fragment of bovine κ-CN. The resulting models were evaluated using the predicted local distance difference test (pLDDT) confidence metric. The pLDDT values ranged from 69.2 to 73.8 (Table S1), indicating moderate confidence in the predicted local structure of this short peptide fragment. Among the generated models, model No. 4 showed the highest pLDDT value and was therefore selected as the starting structure for subsequent calculations.

3.2. Molecular Dynamics Simulations of Chn-SI, Chn-Act, and Chn-Fold

Ten independent MD simulations were performed for each of the three bovine chymosin models: Chn-Act, Chn-SI, and Chn-Fold. The distribution of the Tyr77 torsion angle across the concatenated trajectory frames is shown in Figure 5, while the corresponding RMSD profiles are presented in Figure S1. In all three systems, RMSD values remained below 5 Å.

In the Chn-Fold and Chn-Act simulations, Tyr77 predominantly sampled the open-state torsion-angle region near 300°, which accounted for more than 90% of the analyzed frames. In simulations initiated from the Chn-SI structure, the self-inhibited-state region near 175° was populated for approximately 40% of the analyzed frames, whereas the remaining frames sampled the open-state region near 300°.

Under the present simulation conditions, the Chn-Fold model showed stable RMSD behavior and predominantly sampled the open-state Tyr77 conformation. These features supported its use as the starting model for subsequent conformational analysis and substrate-docking calculations.

For each of the 10 Chn-Fold simulations (MD-1–MD-10), the time-averaged RMSD and the corresponding standard deviation σ(RMSD) were calculated (Figure 6). The MD-10 trajectory showed the lowest mean RMSD (2.123 Å) together with moderate conformational fluctuations (σ = 0.19 Å) and was therefore selected for downstream structural clustering. Clustering was performed using the MD-10 trajectory, and the representative structure of the most populated cluster was subsequently used as the receptor structure for molecular docking.

3.3. Molecular Docking Results

Molecular docking of the bovine κ-CN peptide fragment was performed using the representative Chn-Fold structure obtained from the most populated cluster of the MD-10 trajectory. Docking poses were generated under the positional constraint that the Phe105–Met106 region of the substrate should remain within 7 Å of the Asp34/Asp216 catalytic dyad. Among the resulting solutions (Table S2), a pose showing favorable spatial accommodation of the peptide within the catalytic cleft was selected for subsequent analysis.

In the selected docking pose, the κ-CN peptide occupied the catalytic cleft, with the Phe105–Met106 region oriented toward Asp34 and Asp216 (Figure 7). The distances between the carbonyl oxygen of Phe105 and the carboxyl groups of Asp34 and Asp216 were 2.88 and 4.00 Å, respectively. Because proximity to the catalytic dyad was included as a positional constraint during docking, these distances should be interpreted as geometric characteristics of the selected starting model rather than as independent evidence of a catalytically competent configuration.

Within the selected docking pose, the Phe105 carbonyl group was positioned asymmetrically relative to the catalytic dyad, with distances of 2.88 Å to Asp34 and 4.00 Å to Asp216. This geometry places the scissile Phe105–Met106 region in the vicinity of the catalytic residues. However, because proximity to the Asp34/Asp216 dyad was included as a positional constraint during docking, these distances should not be interpreted as evidence for charge redistribution, transition-state stabilization, or formation of a catalytically competent state.

In the selected Chn-Fold + κ-CN pose, Tyr77 adopted a torsion angle of approximately 300°, corresponding to the open-state geometry and leaving the S1 region accessible to the Phe105 side chain. According to the Schechter–Berger nomenclature [64], Phe105 occupies the P1 position of the substrate and the corresponding enzyme pocket is designated S1. No direct π–π stacking interaction was observed between Tyr77 and Phe105 in this pose. Tyr77 also maintained contacts with neighboring residues, including Trp41, Ile75, and Gly80, which are consistent with retention of the open-state geometry. In addition, a hydrogen-bond interaction between Asp279 of chymosin and His100 of the κ-CN His-Pro cluster was observed (Table S3), which may contribute to positioning of the substrate peptide within the binding cleft.

An analysis of the Ramachandran angles [65,66] for the amino acid residues of both partners of the complex showed that the overwhelming majority of dihedral angles are located in the energetically allowed regions (Figure S2), which indicates the stereochemical correctness of the constructed model. The only mismatch is observed for the Gln280 residue (marked with a green circle in Figure S2), located in the loop region near the binding cleft. This deviation is located in a flexible loop region and does not substantially affect the overall stereochemical quality of the model. The resulting structure demonstrates a high degree of complementarity of the enzyme and substrate surfaces, characteristic of the homologous Chn-Fold + κ-CN system.

Thus, the selected docking pose provides a structurally plausible starting model of the Chn-Fold + κ-CN complex, in which the substrate occupies the catalytic cleft and the Phe105–Met106 region is positioned in the vicinity of the Asp34/Asp216 catalytic dyad. This model was therefore used as the starting structure for subsequent molecular dynamics simulations of the enzyme–substrate complex.

3.4. Molecular Dynamics Simulations of the Chn-Fold + κ-CN Complex

The selected docking pose of the Chn-Fold + κ-CN complex was used as the starting structure for 10 independent MD simulations. The RMSD profiles showed trajectory-dependent fluctuations followed by overall stabilization of the complex (Figure S3).

Analysis of the Tyr77 torsion angle showed that the enzyme predominantly retained the open-state geometry, with values mainly distributed around 280°–300° (Figure 8). No persistent transition to the self-inhibited Tyr77 conformation was observed within the sampled trajectories. These results are consistent with stabilization of the open-state geometry in the substrate-bound complex under the present simulation conditions.

For each of the 10 Chn-Fold + κ-CN simulations (MD-1–MD-10), the time-averaged RMSD and the corresponding standard deviation σ(RMSD) were calculated (Figure 9). The MD-5 trajectory was selected for subsequent structural analysis based on its RMSD behavior and conformational stability. Further analysis of enzyme–substrate interactions and structural clustering was therefore performed using the MD-5 trajectory. The representative structure of the most populated MD-5 cluster was used for subsequent structural characterization of the pre-reactive complex.

3.5. Enzyme–Peptide Interaction

Analysis of enzyme–substrate interactions was performed for the selected 300-ns MD-5 trajectory of the Chn-Fold + κ-CN complex. A persistent network of intermolecular contacts was observed within this trajectory (Figure 10, Table S4), consistent with stabilization of a pre-reactive enzyme–substrate complex. In the N-terminal region of the κ-CN fragment, hydrogen-bond interactions involved His76 (Chn-Fold)–Pro99 (κ-CN), His76 (Chn-Fold)–His102 (κ-CN), Ser132 (Chn-Fold)–His100 (κ-CN), and Asp34 (Chn-Fold)–His102 (κ-CN). These interactions contribute to anchoring of the substrate peptide within the binding cleft. No persistent direct contact between the Phe105 carbonyl group and either Asp34 or Tyr77 was observed in the analyzed MD-5 trajectory, indicating that the scissile region retained local conformational mobility within the bound complex.

Additional contacts were observed in the central and C-terminal regions of the substrate, including hydrogen bonds involving Thr219–Ala107, Ser220–Ile108, Gln280–Pro110, Glu118–Lys111, Thr115–Gln114, Cys52–Lys116, and Phe114–Lys116, together with a Glu118–Lys111 salt bridge (Figure 10, Table S4). Together, these interactions form a multipoint substrate-anchoring network extending across the κ-CN fragment. The absence of a persistent Phe105–Asp34 contact is consistent with a dynamically organized pre-reactive state in which the substrate remains bound while the scissile region retains local mobility relative to the catalytic dyad.

Figure 10. Intermolecular interactions identified in the selected 300-ns MD-5 trajectory of the Chn-Fold + κ-CN complex. The contacts shown correspond to interactions with 100% occupancy within the analyzed MD-5 trajectory. κ-CN residues are shown as molecular structures and chymosin residues as colored ellipsoids according to the legend. Black dotted lines indicate hydrogen bonds, and the red dotted line indicates a salt bridge.

Figure 10. Intermolecular interactions identified in the selected 300-ns MD-5 trajectory of the Chn-Fold + κ-CN complex. The contacts shown correspond to interactions with 100% occupancy within the analyzed MD-5 trajectory. κ-CN residues are shown as molecular structures and chymosin residues as colored ellipsoids according to the legend. Black dotted lines indicate hydrogen bonds, and the red dotted line indicates a salt bridge.

3.6. MD Clustering Analysis

Structural clustering was performed using the selected 300-ns MD-5 trajectory of the Chn-Fold + κ-CN complex. The representative structure of the most populated cluster was selected for subsequent structural analysis and is shown in Figure 11. In this representative structure, several intermolecular hydrogen-bond contacts were observed between chymosin and the κ-CN peptide (Table S5).

In the representative structure of the most populated MD-5 cluster, hydrogen-bond interactions were observed between His76 of chymosin and Pro99 of κ-CN, and between Gly78 of chymosin and His102 of κ-CN. These contacts contribute to stabilization of the N-terminal region of the peptide within the binding cleft. According to the Schechter–Berger nomenclature [64], Phe105 occupies the P1 position of the substrate, whereas Pro101 occupies the P2 position. The His102 and Pro99 residues are part of the His-Pro cluster of κ-CN. In this clustered structure, Phe105 remains conformationally mobile and does not form a persistent direct interaction with Asp34 or Tyr77, indicating that the scissile region is not rigidly fixed within the catalytic site.

On the C-terminal side of the scissile bond, Pro110 of κ-CN, corresponding to the P4′ position, forms a hydrogen-bond interaction with Gln280 of chymosin. This interaction may contribute to stabilization of the C-terminal region of the substrate within the S3′–S4′ region of the binding cleft.

The representative structure of the most populated cluster obtained from the MD-5 trajectory was further analyzed to characterize the relative positioning of the substrate with respect to the Asp34/Asp216 catalytic dyad (Figure 12).

In this representative clustered structure, Tyr77 adopts an open-state conformation with a torsion angle of 304.3°, leaving the catalytic cleft accessible to the substrate peptide. The κ-CN fragment occupies the binding cleft with the Phe105–Met106 region oriented toward the Asp34/Asp216 catalytic dyad.

The distance between the carbonyl oxygen O3 of Phe105 and the O4 atom of Asp34 is 7.05 Å. This distance is substantially larger than a typical hydrogen-bonding distance and indicates that Phe105 is not rigidly engaged with Asp34 in this representative structure. Instead, the scissile region remains spatially associated with the catalytic cleft while retaining local mobility relative to the catalytic dyad, consistent with a dynamically organized pre-reactive configuration.

3.7. QM/MM Refinement of the Catalytic Region

A separate frame from the selected MD-5 trajectory was used as the starting structure for the QM/MM calculations. This frame was chosen based on two geometric criteria relevant to organization of the catalytic region: an open-state Tyr77 torsion angle and a short Phe105–Asp34 distance. The resulting configuration of the Chn-Fold + κ-CN complex is shown in Figure 13a and was used to construct the QM/MM model.

A volumetric occupancy map of water molecules in the vicinity of Phe105, Met106, Asp34, and Asp216 is shown in Figure 13b. A high-occupancy water region was identified approximately 4 Å from the scissile peptide bond. A water molecule located within this region was selected as the initial putative nucleophilic water for construction of the QM/MM model and for the subsequent relaxed scan along the O(water)···C(carbonyl) coordinate.

The QM/MM-refined configuration corresponding to the minimum along the explored O(water)···C(carbonyl) coordinate is shown as Structure I in Figure 14. In this configuration, the C1–N2 distance remains 1.34 Å, consistent with an intact amide bond, whereas the distance between the O4 atom of the selected water molecule and the carbonyl carbon C1 of Phe105 is 2.73 Å. The local geometry is stabilized by a hydrogen-bond network involving the scissile peptide region, Asp34, and nearby water molecules. These structural features are consistent with a pre-reactive configuration of the catalytic region.

Further reduction of the O4···C1 distance along the imposed scan coordinate was accompanied by rearrangement of the local active-site geometry and elongation of the C1–N2 bond. The energetic increase sampled along this one-dimensional relaxed scan reached approximately 12.9 kcal mol⁻¹ relative to the selected minimum (Figure 14). Because no transition-state optimization or vibrational analysis was performed, this value should be interpreted as an approximate energetic penalty associated with progressive water approach along the chosen coordinate rather than as an activation barrier. Accordingly, the higher-energy structure sampled during the scan is not assigned here as a transition state or tetrahedral intermediate.

4. Discussion

In the present study, a multiscale computational approach was used to characterize the structural organization of the bovine Chn/κ-CN complex preceding cleavage of the Phe105–Met106 bond. Particular attention was paid to the conformational behavior of Tyr77, substrate anchoring within the binding cleft, and the local organization of the catalytic region.

The replicated MD simulations revealed substantial conformational plasticity of the Tyr77-containing flap. In simulations initiated from the Chn-Fold and Chn-Act structures, Tyr77 predominantly sampled the open-state torsion-angle region near 300°, accounting for more than 90% of the analyzed frames. In simulations initiated from the Chn-SI structure, the self-inhibited-state region near 175° was populated for approximately 40% of the analyzed frames, whereas the remaining frames sampled the open-state region. These results indicate that both conformational states are accessible under the present simulation conditions, with the open-state geometry being more frequently sampled in the systems investigated here.

This behavior is consistent with previous studies suggesting that chymosin can undergo interconversion between open and self-inhibited conformations [21,28,30,41,67]. At the same time, the present simulations should not be interpreted as a direct determination of the equilibrium populations of these states in solution, because the observed distributions depend on the starting structures, simulation protocol, and accessible timescale. Rather, the results demonstrate that the open Tyr77 configuration is dynamically accessible and can remain stable over the simulated trajectories. More broadly, related studies published in Computation have demonstrated the utility of combining molecular docking and MD simulations for the characterization of biomolecular complexes [68,69], as well as the application of QM/MM approaches to investigate local chemical reactivity in biological systems [70].

In the Chn-Fold + κ-CN simulations, Tyr77 remained predominantly within the open-state angular range, without a persistent transition to the self-inhibited geometry. The bound peptide formed a network of contacts involving residues of the His-Pro region, including interactions with His76, Ser132, Asp34, and other residues of the binding cleft. These observations are consistent with a contribution of substrate contacts to stabilization of the open-state geometry. However, the present simulations do not by themselves establish a unique allosteric mechanism by which the His-Pro cluster drives the Chn-SI-to-Chn-Act transition.

Analysis of the selected 300-ns MD-5 trajectory and of the representative structure of its most populated cluster revealed a multipoint pattern of κ-CN anchoring within the chymosin binding cleft (Figure 10, Figure 11 and Figure 12). Contacts involving the N-terminal His-Pro region included interactions of Pro99 and His102 of κ-CN with residues His76, Gly78, Ser132, and Asp34 of chymosin, whereas the central and C-terminal regions of the substrate formed additional contacts with residues including Glu118 and Gln280. Together, these interactions stabilize several regions of the peptide simultaneously, while allowing the Phe105–Met106 scissile region to retain local conformational mobility.

In the representative structure of the most populated MD-5 cluster, the distance between the carbonyl oxygen of Phe105 and Asp34 was approximately 7.05 Å, indicating the absence of a direct hydrogen-bond interaction between these groups in this configuration (Figure 12). This value differs substantially from the 2.88 Å distance observed in the selected docking pose. However, because proximity of the Phe105–Met106 region to the Asp34/Asp216 dyad was included as a positional constraint during docking, this difference should not be interpreted simply as disruption of a pre-existing catalytic contact. Rather, the MD results indicate that the scissile region can move relative to the catalytic dyad while the peptide remains anchored through contacts elsewhere along the binding cleft.

In the selected docking configuration shown in Figure 7, the aromatic rings of Tyr77 and Phe105 were spatially separated and did not adopt a geometry compatible with direct π–π stacking. Thus, in this particular open-state configuration, Tyr77 does not directly anchor the Phe105 side chain. Instead, the open orientation of Tyr77 appears to maintain accessibility of the catalytic cleft while other enzyme–substrate contacts contribute to peptide positioning.

Taken together, the MD results support a model in which κ-CN binding is stabilized by a distributed network of contacts rather than by rigid fixation of the scissile bond itself. Such multipoint anchoring may help maintain the substrate within the binding cleft while preserving local flexibility of the Phe105–Met106 region, thereby providing a dynamically organized pre-reactive ensemble from which geometries more favorable for catalytic organization can be sampled.

QM/MM refinement of the catalytic region provided a more detailed description of the local organization of the selected pre-reactive configuration. The optimized structure revealed a hydrogen-bonded network involving Asp34, Asp216, the Phe105–Met106 region, and several nearby water molecules. Rather than identifying a single rigidly preorganized catalytic water molecule, the calculations indicate that solvent molecules form part of a dynamically organized hydrogen-bond network surrounding the scissile bond. Such an arrangement is compatible with proposed mechanisms of aspartic proteases in which catalytic organization involves concerted rearrangement of the substrate, catalytic aspartates, and active-site water molecules [37,71,72,73].

The exploratory one-dimensional QM/MM relaxed scan further showed that progressive approach of the selected water molecule toward the carbonyl carbon of Phe105 is energetically accessible within the modeled active-site environment. The energetic increase along the imposed O(water)···C(carbonyl) coordinate reached approximately 12.9 kcal mol⁻¹ relative to the selected minimum. Because the calculation was not accompanied by transition-state optimization or vibrational analysis, this value should not be interpreted as an activation barrier, and the higher-energy geometries sampled along the scan cannot be assigned unambiguously to a transition state or tetrahedral intermediate. Instead, the scan provides a qualitative measure of the energetic response of the active site to nucleophilic water approach and supports the structural plausibility of further catalytic rearrangement from the modeled pre-reactive state.

5. Conclusions

This study combined AlphaFold3-based structural modeling, molecular docking, replicated molecular dynamics simulations, and QM/MM calculations to characterize the pre-reactive complex formed by bovine chymosin and the Chn-sensitive fragment of κ-casein. Replicated MD simulations showed that the open-state Tyr77 geometry is dynamically accessible under the simulated conditions and remains predominant in the Chn-Fold + κ-CN complex. Analysis of the selected MD-5 trajectory revealed a multipoint substrate-anchoring network involving both the N- and C-terminal regions of the κ-CN fragment, while the Phe105–Met106 scissile region retained local mobility relative to the Asp34/Asp216 catalytic dyad.

QM/MM refinement of a selected MD-5 configuration revealed an organized hydrogen-bond network involving Asp34, Asp216, the Phe105–Met106 region, and nearby water molecules. An exploratory one-dimensional relaxed scan along the O(water)···C(carbonyl) coordinate showed that progressive approach of the selected water molecule is associated with a moderate energetic increase of approximately 12.9 kcal mol⁻¹ relative to the selected minimum. Because no transition-state optimization or vibrational analysis was performed, this value should not be interpreted as an activation barrier. Instead, the scan supports the structural and energetic plausibility of nucleophilic water approach within the modeled active-site environment.

Overall, the results support a model in which bovine κ-CN recognition is achieved through distributed enzyme–substrate contacts, whereas the scissile Phe105–Met106 region remains dynamically organized rather than rigidly fixed before catalysis. The resulting multiscale model provides an atomistic framework for further investigation of chymosin substrate specificity and catalytic organization. A complete description of the hydrolysis pathway, including transition-state characterization and determination of the corresponding activation free energy, will require dedicated reaction-path calculations beyond the scope of the present study.