The X-ray crystal structure of DAPK1 that contains the CD and ARD domains in complex with Ca²⁺/CaM was obtained from the Protein Data Bank (PDB ID: 2X0G) (de Diego et al., 2010). The coordinates of the missing residues of CaM (residues 74–84) at the flexible linker were built using the MODELLER program (Fiser and Sali, 2003). All Ca²⁺ were removed from CaM to simulate the Ca²⁺-unbound DAPK1−CaM system. The protonation states for the ionizable residues at pH = 7 were used. For the histidine residues, the protonation states were determined based on the PROPKA calculation (Rostkowski et al., 2011). All the histidine residues in the DAPK1−CaM complex were set in a neutral state, however, the HID or HIE forms for the histidine residue was chosen based on the local hydrogen bonding network.

The tleap module of AMBER 18 was used to add the hydrogen atoms for the DAPK1−CaM complex (Salomon-Ferrer et al., 2013). The force field parameters for Ca²⁺ were downloaded from the AMBER parameters database (http://www.pharmacy.manchester.ac.uk/bryce/amber), using the value of van der Waals radius of R^* = 1.79 Å and the well depth of ε = 0.0140 kcal/mol. This force field parameter for the Ca²⁺ used in the MD simulations have been reported to reproduce the Ca²⁺−residue interactions in the X-ray crystal structure (Lawrenz et al., 2010).

Both systems were embedded in a truncated octahedron box of the TIP3P water molecules with a 10 Å buffer (Jorgensen et al., 1983). The Na⁺ counter-ions were added to maintain the electroneutrality of both systems, and then 0.15 mol/L NaCl was added to simulate physiological conditions. Energy minimizations and MD simulations were performed in an isothermal isobaric ensemble (NPT) with periodic boundary conditions using the CUDA module. Structure optimizations were carried out in a stepwise manner with the harmonic force restraint of the DAPK1−CaM complex and Ca²⁺ and the whole system then allowed to fully move. Energy optimizations were performed using the steepest descent method for the first 5,000 steps and then the conjugated gradient method for the subsequent 10,000 steps. Then, both systems were heated from 0 to 300 K in 100 ps. This was followed by constant temperature equilibration at 300 K for 500 ps at the canonical ensemble (NVT). In the production runs, 1,000 ns simulations were carried out, which were repeated three times using random velocities. The long-range electrostatic interactions were treated by the particle mesh Ewald method (Darden et al., 1993). A cut-off of 10 Å was used for short-range electrostatics and van der Waals interactions. All covalent bonds involving hydrogen atoms were constrained using the SHAKE method (Ryckaert et al., 1977). An integration step of 2 fs was used.

The cross-correlation matrix (C _ij) between the fluctuations of the Cα atoms of the DAPK1–CaM was used to show the coupling of the motions between the protein residues (Li et al., 2020, Li et al., 2021; Shi et al., 2022). C _ij was calculated using following equation, C i , j = c i , j c i , i 1 / 2 c j , j 1 / 2

Positive C _ij values mean the two atoms i and j moving in the same direction, while negative C _ij values describe anti-correlated motions between the two atoms i and j.

The correlation C _ij data were further used to weight edges and compute the edge distance D i j using the following equation (Sethi et al., 2009; Saha et al., 2020; Zhuang et al., 2022a; 2022b), which represents the possibility of information flow. D i j = − log C i j

The community network was defined as a set of nodes (Cα atoms of the DAPK1–CaM complex), which was described by CC _ij weighted edges between two nodes within a cutoff distance of 4.5 Å for >75% of the simulated trajectory. The distribution of communities within the whole protein was decided and optimized using the Girvan-Newman algorithm. Communities harboring residues less than three were omitted.

The molecular mechanisms–generalized Born surface area energy calculations (MM–GBSA) were performed using the following equations (Wang et al., 2019; Liu et al., 2022a; Liu et al., 2022b; He et al., 2022; Zhu et al., 2022). ∆ G b i n d i n g = ∆ G c o m p l e x − ∆ G p r o t e i n + ∆ G l i g a n d ∆ G b i n d i n g = ∆ G E g a s + ∆ G s o l v a t i o n − T ∆ S ∆ E g a s = ∆ E v d W + ∆ E e l e ∆ G s o l v a t i o n = ∆ G G B + ∆ G n o n p o l a r ∆ G n o n p o l a r = γ × S Α S Α + b

The ΔE _gas, ΔE _vdW, ΔE _ele, ΔG _solvation, ΔG _GB, and ΔG _nonpolar terms were gas energy, van der Waals energy, electrostatic energy, solvation free energy, the polar energy, and the non-polar energy, respectively. The ΔG _nonpolar was calculated using the function of the solvent accessible surface area (SASA) with the γ value of 0.0072 kcal/(mol Ų) and the b value of 0 kcal/mol. The TΔS energy item was not calculated due to the extremely long durations of normal mode analysis for large systems.

To investigate the overall system stability and time-dependent conformational dynamics of the DAPK1−CaM complex in the presence and absence of Ca²⁺, a total of 6 μs-length MD simulations were carried out in the explicit water environment using the AMBER 18 program. Both systems, including Ca²⁺-bound and-unbound DAPK1–CaM complexes, were simulated in three independent times using random velocities to achieve reliable statistics. The overall stability of protein complex and the time-dependent protein-protein interactions were analyzed through extracting structural complexes from the production trajectories at different time intervals (Lu et al., 2019, 2021; Jang et al., 2020; Maloney et al., 2021; Ni et al., 2021). Then, conformational changes of the DAPK1–CaM complex in the Ca²⁺-bound and-unbound forms were compared. Root-mean-square deviation (RMSD), principal component analysis (PCA), dynamical cross-correlation matrix (DCCM), radial pair distribution function g(r), community network analysis, and MM–GBSA binding free energy calculations were calculated throughout 1,000 ns time scale simulations.

To reveal the system convergence, the RMSD of the Cα atoms for the complex, DAPK1, and CaM was monitored with reference to the crystal structure (PDB ID: 2X0G) along the three independent simulations. In general, the complex (Figure 2A), DAPK1 (Figure 2B), and CaM (Figure 2C) in both systems were convergent in the early stage of simulations. Moreover, the RMSD values for the complex, DAPK1, and CaM were similar in both systems. For instance, in the Ca²⁺-bound state, the RMSD values for the complex, DAPK1, and CaM were 3.29 ± 0.39, 1.49 ± 0.18, and 4.04 ± 0.42 Å, respectively. In the Ca²⁺-unbound state, the RMSD values for the complex, DAPK1, and CaM were 3.02 ± 0.30, 1.51 ± 0.16, and 3.41 ± 0.41 Å, respectively. Together, these results suggested that the system stability of the DAPK1–CaM had no significant change in the absence of Ca²⁺ binding.

To reveal the large-scale collective motions and the conformational interconversion of the DAPK1−CaM complex in the Ca²⁺-bound and-unbound states, PCA analysis of both simulated systems was performed (Masterson et al., 2011; An et al., 2021). The covariance matrix of the Cα atoms for the DAPK1−CaM complex was diagonalized to produce a set of eigenvalues and the corresponding eigenvectors. Each eigenvector is named the principal component (PC), which is related to an eigenvalue corresponding to the mean square fluctuation projected along the that eigenvector. The first several PCs represent the overall fluctuations of the structural complex (Palermo et al., 2016). Each snapshot was first subjected to RMS-fit to the initial crystal structure of the DAPK1−CaM complex (PDB ID: 2X0G) as the same reference configuration. Then, all snapshots from MD simulations were projected into the collective coordinate space defined by the first two eigenvectors (PC1 and PC2), which reflected the essential conformational subspace sampled by the DAPK1−CaM complex in the Ca²⁺-bound and-unbound states.

We performed PCA analysis of the DAPK1−CaM complex in the presence and absence of Ca²⁺ binding and found that the PC1 and PC2 represented ∼60% of variance in coordinates along the MD simulations in both systems. The free energy landscapes (FELs) of the PC1 and PC2 showed the distinct conformational space sampled by the DAPK1−CaM complex in the Ca²⁺-bound and-unbound states. As shown in Figure 3A, in the presence of Ca²⁺, the PC1 and PC2 plots sampled a broad distribution in the FELs, with the PC1 and PC2 values in the range of ∼ − 60 to ∼60 and ∼ − 40 to ∼60, respectively. However, in the absence of Ca²⁺ (Figure 3B), the PC1 and PC2 plots sampled a confined distribution, with the PC1 and PC2 values in the range of ∼ − 30 to ∼40 and ∼ − 30 to ∼20, respectively. Moreover, two major conformational substates of the DAPK1−CaM complex were observed in both the Ca²⁺-bound and-unbound states. Collectively, these results suggested that Ca²⁺ binding increased the dynamics and the conformational space of the DAPK1−CaM complex.

We next sought to explore how Ca²⁺ induces conformational changes of the DAPK1−CaM complex. The global intra- and inter-domain motions in the protein-protein complex involve the collective motions of the protein backbone atoms. To reveal the correlated or anti-correlated motions between different residues and domains within the DAPK1−CaM complex, the dynamical cross-correlation matrix (DCCM) for the Cα atoms of the structural complex were constructed and analyzed in both the Ca²⁺-bound and-unbound states. As shown in Figure 4, the DCCM plots showed the correlated or anti-correlated motions between residues within the protein complex, with the regions in orange representing the correlated motions and the regions in purple representing the anti-correlated motions and the density of the color corresponding to the intensity of correlated/anti-correlated motions.

In the presence of Ca²⁺ (Figure 4A), the intra-domain motions in both the independent DAPK1 and CaM showed mainly correlated motions, while the inter-domain motions between the DAPK1 and CaM underwent weak anti-correlated motions. However, in the absence of Ca²⁺ (Figure 4B), the intra-domain motions in the independent DAPK1 showed mixed correlated and anti-correlated motions. Although the intra-domain motions in the independent DAPK1 showed correlated motions without Ca²⁺ binding, the intensity of the correlated motions was stronger in the Ca²⁺-unbound state than in the Ca²⁺-bound state. Remarkably, compared to the inter-domain motions between the DAPK1 and CaM in the Ca²⁺-bound state, they underwent enhanced anti-correlation motions in the Ca²⁺-unbound state. The increase of the anti-correlation motions between the DAPK1 and CaM in the Ca²⁺-unbound state indicated that the Ca²⁺-unbound CaM might disassociate from the DAPK1, which can be assessed by the following binding free energy calculations.

To explore the effect of Ca²⁺ removal on the allosteric network of the DAPK1–CaM complex, the difference contact network analysis (dCNA) was performed (Li et al., 2022). The connected residues are considered as one community, which can be served as a synergistic functional unit of the protein. The dCNA analysis can unravel the intercommunity changes between difference conformational ensembles. The same communities are represented by spheres that are connected by sticks, whose width is corresponding to the intensity of interaction between different communities. Figure 5 shows the three-dimensional and two-dimensional visualizations of the dCNA between the Ca²⁺-bound and-bound states, respectively. We mainly focused on the connectivity between the interface of the DAPK1–CaM complex. The community one mainly contained the B-loop and the P-loop of the DAPK1, the community 3 was composed of the C-terminal ARD of the DAPK1, the community six included the helices α1 and α4 of the CaM(N), and the community 8 possessed the CaM(C) that contained the helices α5–α8. After the removal of Ca²⁺, the connectivity between the community 1 of the DAPK1 and the community six of the CaM(N), between the community 3 of the DAPK1 and the communities 6/7 of the CaM(N), and the community 4 of the DAPK1 and the community 7 of the CaM(N) decreased. These results suggested that the allosteric communication between the DAPK1 and CaM was reduced in the Ca²⁺-bound state, which would weaken the DAPK1–CaM interaction.

To further show the energetics of the DAPK1–CaM interactions in the Ca²⁺-bound and-unbound states, the binding free energy (ΔG_binding) calculations were performed using the MM–GBSA method, which has been proved successfully to evaluate protein-ligand or protein-protein interactions (Liu et al., 2018; Wang et al., 2019, 2021, 2022). A total of 200 snapshots for the last 200 ns snapshots were selected from the MM–GBSA binding free energy calculations. As shown in Table 1, the ΔG_binding values in the Ca²⁺-bound and-unbound states were −112.98 ± 20.63 and −86.40 ± 18.79 kcal/mol, respectively. Thus, the binding free energy between DAPK1 and CaM was larger by −26.58 kcal/mol in the Ca²⁺-bound state than that in the Ca²⁺-unbound state, which suggested that the removal of Ca2+ would disfavor the DAPK1–CaM interactions.

Ca²⁺ has a radius of 0.99 Å, showing coordination flexibility with coordination numbers (CNs) ranging from six to eight and with an average CN of 7.3 (Kirberger et al., 2008). In the X-ray crystal structure of the DAPK1−CaM complex (PDB ID: 2X0G) (de Diego et al., 2010), there are four Ca²⁺ binding sites, including two in the C-terminus (CS1 and CS2) and two in the N-terminus (NS1 and NS2). However, during crystallization, crystal waters are often omitted. With this idea in mind, we calculated the radial pair distribution function g(r) between each Ca²⁺ and the water molecules. Radial pair distribution function g(r) can provide the probability of finding particles at a certain distance (Lu et al., 2013; Sameera et al., 2020). Figure 6 shows the radial pair distribution function g(r) derived from the MD simulations in the Ca²⁺-bound DAPK1−CaM complex. For the Ca²⁺ binding sites at the CS1, CS2, and NS1, all plots showed a sharp peak at 2.35 Å, which suggested that there was a high probability that the water molecules were coordinated to these Ca²⁺ during MD simulations. However, in the NS2 binding site, no sharp peak at 2.35 Å was observed, indicating that no water molecules were coordinated to Ca²⁺ in the NS2 binding site. Furthermore, the CNs of the water molecules for Ca²⁺ at the CS1, CS2, and NS1 binding sites were computed. As shown in Figure 6, one water molecule was coordinated to each Ca²⁺ at the CS1, CS2, and NS1 binding sites.

To further characterize the Ca²⁺ coordination mode in each binding site, a detailed analysis of the MD trajectories was performed using cluster analysis. Figure 7 shows the most representative structure obtained from the MD trajectories and the detailed coordination mode for each Ca²⁺ binding site. Each Ca²⁺ has a CN of 7 from the CaM. In the CS1 binding site (Figure 7A), Ca²⁺ was coordinated by one oxygen atom from the sidechain Asn97, one oxygen atom from the sidechain Asp93, two oxygen atoms from the sidechain Asp95, two oxygen atoms from the sidechain Glu104, and one water oxygen atom. In the CS2 binding site (Figure 7B), Ca²⁺ was coordinated by one oxygen atom from the sidechain Asp129, one oxygen atom from the sidechain Asp131, one oxygen atom from the sidechain Asp133, one oxygen atom from the sidechain Gln135, two oxygen atoms from the sidechain Glu140, and one water oxygen atom. In the NS1 binding site (Figure 7C), Ca²⁺ was coordinated by one oxygen atom from the sidechain Asp20, one oxygen atom from the sidechain Asp22, one oxygen atom from the sidechain Asp24, one oxygen atom from the backbone Thr26, two oxygen atoms from the sidechain Glu31, and one water oxygen atom. In the NS2 binding site (Figure 7D), Ca²⁺ was coordinated by one oxygen atom from the sidechain Asp56, one oxygen atom from the sidechain Asp58, one oxygen atom from the sidechain Asn60, one oxygen atom from the backbone Thr62, one oxygen atom from the sidechain Asp64, and two oxygen atoms from the sidechain Glu67. These exquisite coordinated networks between Ca²⁺ and CaM render the CaM in a restrained conformational change when it binds to the DAPK1 protein. Thus, CaM would undergo large conformational changes after the removal of Ca²⁺ from each binding site, which was unfavorable for the binding to DAPK1 as revealed by the MM−GBSA binding free energy calculations.

To investigate how Ca²⁺ removal affected the detailed conformational changes of the CaM and the DAPK1−CaM complex, the most representative structural complexes in both the Ca²⁺-bound and-unbound states were extracted using cluster analysis of MD trajectories (Shao et al., 2007). We first overlapped the independent CaM from MD simulations in both systems. The RMSD of the CaM between the Ca²⁺-bound and-unbound states was 2.07 Å. As shown in Figure 8, compared to the Ca²⁺-bound state, the appreciable conformational rearrangements of the CaM were found in the helices α1, α3, α4, α7 and the loop that connects the helices α1 to α2. According to the X-ray crystal structure of the DAPK1−CaM complex (PDB ID: 2X0G), the helices α1, α4, α7 and the loop that connects the helices α1 to α2 from the CaM are involved in the interactions with both the CD and the ARD of DAPK1. Thus, the disturbed conformational changes at the CaM interface would allosterically affect the association with the DAPK1. As a result, Ca²⁺ may serves an allosteric modulator of DAPK1−CaM protein-protein interaction (PPI), which has been an area of intensive research for PPI drug discovery (Ni et al., 2019, 2022).

To further reveal the effect of Ca²⁺ removal on the conformational dynamics of the DAPK1−CaM complex, we superimposed the structural complexes in both the Ca²⁺-bound and-unbound states. As shown in Figure 9, compared to the DAPK1−CaM in the Ca²⁺-bound state, the removal of Ca²⁺ allosterically changed the conformations of the DAPK1 at the basic loop (B-loop) and the phosphate-binding loop (P-loop), while the C-terminal ARD of the DAPK1 underwent little conformational changes. Indeed, the B-loop and the P-loop of the DAPK1 are involved in the interactions with the helices α1 and α4 and the loop connecting the helices α1 to α2 of the CaM. The C-terminal ARD of the DAPK1 engages with the helices α2 and α7. Therefore, the disruption of the interface of the DAPK1-CaM complex after Ca²⁺ removal would disassociate the binding of Ca²⁺-free CaM to the DAPK1, which has been revealed by the above binding free energy calculations. As a result, the catalytic efficiency of DAPK1 would reduce without the regulation of CaM.