SMD Simulations Research Articles

A comparative study on penetration mechanisms of drug-loaded carbon and boron nitride nanotubes through biological membranes by steered molecular dynamics simulations

Understanding the mechanism by which nanotubes penetrate cell membranes is challenging from multiple perspectives. As a drug delivery system, boron nitride nanotubes (BNNTs) have a similar structure to carbon nanotubes but with B–N bonds instead of C–C bonds. Through a remarkable series of direct and indirect observations within specific cells, these nanotubes are popularly attributed as superior penetrant into cell membrane. Suitable functional groups and polymers are often needed to enhance the biocompatibility and solubility of BNNTs and CNTs in biological media. In addition, to figure out the effect of functional groups, the nanostructures without functional groups were first examined together with their anticancer drugs (fluorouracil and letrozole). All partial charges of the drug and nanotube have been investigated through population analysis. After that, with a total of 40 simulations (MD and SMD simulations), various analytical techniques were employed to examine the interaction between drugs and nanotubes with POPE, which is a class of phospholipids existing in biological membrane, in aqueous media. Noteworthy among these techniques is the mean-squared displacement analysis to compare the diffusion rate of nanocarriers and the radial distribution function analysis, which was utilized to compare water concentrations surrounding nanotubes. Additionally, the stability of the drug within the nanotube was assessed through mass center distance analysis. The diffusion coefficients of the nanotube-membrane complex were compared against various chemical agents by employing mean squared displacement analysis. The findings of the study revealed that the tethering of tetra ethylene glycol results in the augmentation of the water molecules surrounding the nanotubes while simultaneously enhancing the durability of the drug being conveyed. Nonetheless, the addition of tetra ethylene glycol resulted in a reduction in the nanocarrier’s diffusion coefficient. Communicated by Ramaswamy H. Sarma

Different aspects in explaining how mutations could affect the binding mechanism of receptor binding domain of SARS-CoV-2 spike protein in interaction with ACE2.

During replication, some mutations occur in SARS-CoV-2, the causal agent of COVID-19, leading to the emergence of different variants of the virus. The mutations that accrue in different variants of the virus, influence the virus' ability to bind to human cell receptors and ability to evade the human immune system, the rate of viral transmission, and effectiveness of vaccines. Some of these mutations occur in the receptor binding domain (RBD) of the spike protein that may change the affinity of the virus for the ACE2 receptor. In this study, several in silico techniques, such as MD and SMD simulations, were used to perform comparative studies to deeply understand the effect of mutation on structural and functional details of the interaction of the spike glycoprotein of SARS-CoV-2, with the ACE2 receptor. According to our results, the mutation in the RBD associated with the Omicron variant increase binding affinity of the virus to ACE2 when compared to wild type and Delta variants. We also observed that the flexibility of the spike protein of the Omicron variant was lower than in comparison to other variants. In summary, different mutations in variants of the virus can have an effect on the binding mechanism of the receptor binding domain of the virus with ACE2.

Deconstructing allostery by computational assessment of the binding determinants of allosteric PTP1B modulators

Fragment-based drug discovery is an established methodology for finding hit molecules that can be elaborated into lead compounds. However it is currently challenging to predict whether fragment hits that do not bind to an orthosteric site could be elaborated into allosteric modulators, as in these cases binding does not necessarily translate into a functional effect. We propose a workflow using Markov State Models (MSMs) with steered molecular dynamics (sMD) to assess the allosteric potential of known binders. sMD simulations are employed to sample protein conformational space inaccessible to routine equilibrium MD timescales. Protein conformations sampled by sMD provide starting points for seeded MD simulations, which are combined into MSMs. The methodology is demonstrated on a dataset of protein tyrosine phosphatase 1B ligands. Experimentally confirmed allosteric inhibitors are correctly classified as inhibitors, whereas the deconstructed analogues show reduced inhibitory activity. Analysis of the MSMs provide insights into preferred protein-ligand arrangements that correlate with functional outcomes. The present methodology may find applications for progressing fragments towards lead molecules in FBDD campaigns.

Computational drug repurposing effort for identifying novel hits against Y box binding protein 1 as a targeted therapy for ovarian cancer.

Y box binding protein 1 (YB-1) which regulates DNA and RNA associated events is a crucial oncogenic target and its amplified levels have been found in many cancers. YB-1 increases the stability of short-lived mRNAs for multiple oncogenic proteins promoting cancer cell stemness and facilitating the cell cycle's progression at the G1/S phase. The scientific reports of the past three decades strongly suggests that YB-1 may be a potential vital target for treating many drug resistance cancers including ovarian cancer. Thus, in the current study, we conducted a virtual drug repurposing study. For this aim, different computational models were constructed: (i) target-protein driven models (PDB: 6KUG); (ii) machine learning (ML)-based quantitative structure-activity relationships (QSAR) models; (iii) target-protein driven pharmacophore (E-pharmacophore) models. Constructed and statistically validated models were then used to screen FDA approved drugs and compounds in clinical investigation libraries which includes around 10.000 molecules. Successful candidates were used in docking and all-atom molecular dynamics (MD) simulations to further characterize their binding events at the active site of the YB-1 and compute average binding free energies using MM/GBSA. Steered MD (sMD) simulations were also conducted for final hits as well as positive control (SU056). To estimate the binding strength of selected hit compounds to YB-1, sMD simulations were conducted to pull the ligands away from protein at a constant velocity of 10 nm/ns along the vector pointing away from the binding interface. Subsequently, the PMF curves could be used to extract the binding free energies. In this study, the PMF curves were generated by using independent sampling windows from sMD trajectories. Our discovered potential small molecule YB-1 inhibitors can be considered as potential novel inhibitors that can target YB-1.

Cognate RNA-Binding Modes by the Alternative-Splicing Regulator MBNL1 Inferred from Molecular Dynamics.

The muscleblind-like protein family (MBNL) plays a prominent role in the regulation of alternative splicing. Consequently, the loss of MBNL function resulting from sequestration by RNA hairpins triggers the development of a neuromuscular disease called myotonic dystrophy (DM). Despite the sequence and structural similarities between the four zinc-finger domains that form MBNL1, recent studies have revealed that the four binding domains have differentiated splicing activity. The dynamic behaviors of MBNL1 ZnFs were simulated using conventional molecular dynamics (cMD) and steered molecular dynamics (sMD) simulations of a structural model of MBNL1 protein to provide insights into the binding selectivity of the four zinc-finger (ZnF) domains toward the GpC steps in YGCY RNA sequence. In accordance with previous studies, our results suggest that both global and local residue fluctuations on each domain have great impacts on triggering alternative splicing, indicating that local motions in RNA-binding domains could modulate their affinity and specificity. In addition, all four ZnF domains provide a distinct RNA-binding environment in terms of structural sampling and mobility that may be involved in the differentiated MBNL1 splicing events reported in the literature.

Tension Enhances the Binding Affinity of β1 Integrin by Clamping Talin Tightly: An Insight from Steered Molecular Dynamics Simulations.

Integrin activation is a predominant step for cell-cell and cell-ECM interactions. Talin and Kindlin are mechanosensitive adaptor proteins that bind to the integrin cytoplasmic tail and mediate integrin activation, cytoskeleton rearrangement, and focal adhesion assembly. However, knowledge about how Talin and Kindlin synergistically assist integrin activation remains unclear. Here, we performed so-called "ramp-clamp" SMD simulations, which modeled the mechanosignaling from Kindlin, to investigate the effect of tension on the interaction of the β1 integrin cytoplasmic tail with the Talin-F3 domain. The present results showed that mild but not excessive stretching enhanced the binding of integrin with Talin. This mechanical regulation on integrin affinity to Talin referred to an event cascade, in which under stretching, the integrin cytoplasmic tail adopted allostery in response to the mechanical stimulus, remodeling of integrin in favor of Talin-association ensued, and finally, a stable, close-knit complex was formed. In the cascade, the torsion angle transition of integrin was the cue for the stable interaction of the complex under tensile force. The present work suggested a model for Talin and Kindlin to synergistically activate integrin. It should help understand integrin activation and its mechanochemical regulation mechanism, integrin-related innate cellular immune responses, cell adhesion, cell-cell interaction, and integrin-related drug development.

Unveiling molecular interactions that stabilize bacterial adhesion pili

Adhesion pili assembled by the chaperone-usher pathway are superelastic helical filaments on the surface of bacteria, optimized for attachment to target cells. Here, we investigate the biophysical function and structural interactions that stabilize P pili from uropathogenic bacteria. Using optical tweezers, we measure P pilus subunit-subunit interaction dynamics and show that pilus compliance is contour-length dependent. Atomic details of subunit-subunit interactions of pili under tension are shown using steered molecular dynamics (sMD) simulations. sMD results also indicate that the N-terminal “staple” region of P pili, which provides interactions with pilins that are four and five subunits away, significantly stabilizes the helical filament structure. These data are consistent with previous structural data, and suggest that more layer-to-layer interactions could compensate for the lack of a staple in type 1 pili. This study informs our understanding of essential structural and dynamic features of adhesion pili, supporting the hypothesis that the function of pili is critically dependent on their structure and biophysical properties.

Molecular dynamics simulation of calcium silicate hydrate/tannic acid interfacial interactions at different temperatures: configuration, structure and dynamic

In the paper, the effect of temperature on tannic acid molecular conformation and the interaction between tannic acid and substrate were investigated systematically by molecular dynamics. The tensile force of tannin acid separation from CSH matrix at three temperatures was calculated by SMD simulation, and then the adhesion energy between tannin acid and CSH was calculated by Bell Model. The results showed that high temperature affects the gyration of the single bond in the tannic acid, which leads to the increase of tannic acid molecular chain size and radius of gyration. Tannic acid molecules are more extended and the contact area with the substrate became larger in high temperature, which enhanced interfacial adhesion energy. Further, by analyzing the structure and dynamics of the composite, it comes out that not only the substrate has more sites available for tannic acid but also the interfacial ionic bonding formed is more stable after high temperature. Enhanced interfacial chemical bonds limit the bending motion of tannic acid molecular, which leave them in a well-stretched state, to causes tannic acid to be stably adsorbed on the substrate.

Molecular Mechanism for the Allosteric Inhibition of the Human Serotonin Transporter by Antidepressant Escitalopram.

Human serotine transporter (hSERT) is one of the most influential drug targets, and its allosteric modulators (e.g., escitalopram) have emerged to be the next-generation medication for psychiatric disorders. However, the molecular mechanism underlying the allosteric modulation of hSERT is still elusive. Here, the simulation strategies of conventional (cMD) and steered (SMD) molecular dynamics were applied to investigate this molecular mechanism from distinct perspectives. First, cMD simulations revealed that escitalopram's binding to hSERT's allosteric site simultaneously enhanced its binding to the orthosteric site. Then, SMD simulation identified that the occupation of hSERT's allosteric site by escitalopram could also block its dissociation from the orthosteric site. Finally, by comparing the simulated structures of two hSERT-escitalopram complexes with and without allosteric modulation, a new conformational coupling between an extracellular (Arg104-Glu494) and an intracellular (Lys490-Glu494) salt bridge was identified. In summary, this study explored the mechanism underlying the allosteric modulation of hSERT by collectively applying two MD simulation strategies, which could facilitate our understanding of the allosteric modulations of not only hSERT but also other clinically important therapeutic targets.

Dissociation Mechanisms of G-actin Subunits Govern Deformation Response of Actin Filament.

Actin molecules are essential structural components of the cellular cytoskeleton. Here, we report a comprehensive analysis of F-actin's deformation behavior and highlight underlying mechanisms using steered molecular dynamics simulations (SMD). The investigation of F-actin was done under tension, compression, bending, and torsion. We report that the dissociation pattern of conformational locks at intrastrand and interstrand G-actin interfaces regulates the deformation response of F-actin. The conformational locks at the G-actin interfaces are portrayed by a spheroidal joint, interlocking serrated plates' analogy. Further, the SMD simulation approach was utilized to evaluate Young's modulus, flexural rigidity, persistent length, and torsional rigidity of F-actin, and the values obtained were found to be consistent with available experimental data. The evaluation of the mechanical properties of actin and the insight into the fundamental mechanisms contributing to its resilience described here are necessary for developing accurate models of eukaryotic cells and for assessing cellular viability and mobility.

An Investigation of the Influenza Hemagglutinin Membrane Fusion Process using Microsecond-Level MD Simulations

Influenza virus hemagglutinin (HA) is a homotrimeric envelope protein that undergoes pH-triggered conformational changes in order to facilitate the fusion of the viral envelope with the endosomal membrane of the infected host cell. Previous experimental studies have outlined a mechanism for the conformational transitions that lead to membrane fusion. At low pH, the HA1 chains move away from the HA2 chains, which then undergo a series of conformational changes to expose the N-terminal fusion peptide. The HA2 chains then bend at a hinge region to form a hairpin shaped structure, thus bringing the viral envelope into contact with the endosomal membrane. While this model is well established, the conformational changes have never been observed in all-atom equilibrium MD simulations. Using a combination of equilibrium MD simulations and non-equilibrium SMD simulations, we have characterized the conformational dynamics of the membrane fusion process. Six models were simulated for 2.4 microseconds each. By modifying the protonation state of a conserved histidine in the HA2 hinge region, we have successfully simulated the pH-mediated conformational changes that occur when HA1 moves away from HA2. Our atomistic MD simulations confirmed that the HA1 chain associated with an HA2 chain moves away only when the histidine is protonated. Using non-equilibrium pulling simulations, we have also modeled the complex conformational and secondary-structural transitions that lead to the uncoiling of HA2 and the exposure of the fusion peptide. In addition to the conformational changes that occur within hemagglutinin, we have independently modeled the fusion of the viral envelope with the host cell endosomal membrane. For the very first time, we have studied the conformational dynamics of the entire HA-membrane complex using microsecond-level atomistic MD simulations and non-equilibrium simulations.

Enhanced Sampling Simulations of Ligand Unbinding Kinetics Controlled by Protein Conformational Changes.

Understanding unbinding kinetics of protein-ligand systems is of great importance for the design of ligands with desired specificity and safety. In recent years, enhanced sampling techniques have emerged as effective tools for studying unbinding kinetics of protein-ligand systems at the atomistic level. However, in many protein-ligand systems, the ligand unbinding processes are strongly coupled to protein conformational changes and the disclosure of the hidden degrees of freedom closely related to the protein conformational changes so that sampling is enhanced over these degrees of freedom remains a great challenge. Here, we show how potential-scaled molecular dynamics (sMD) and infrequent metadynamics (InMetaD) simulation techniques can be combined to successfully reveal the unbinding mechanism of 3-(1,4-diazabicyclo[3.2.2]nonan-4-yl)-6-[18F]fluo-rodibenzo[b,d]thiophene 5,5-dioxide ([18F]ASEM) from a chimera structure of the α7-nicotinic acetylcholine receptor. By using sMD simulations, we disclosed that the "close" to "open" conformational change of loop C plays a key role in the ASEM unbinding process. By carrying out InMetaD simulations with this conformational change taken into account as an additional collective variable, we further captured the key states in the unbinding process and clarified the unbinding mechanism of ASEM from the protein. Our work indicates that combining sMD and InMetaD simulation techniques can be an effective approach for revealing the unbinding mechanism of a protein-ligand system where protein conformational changes control the unbinding process.

Toward understanding calmodulin plasticity by molecular dynamics.

Aim: Calmodulin interacts in many different ways with its ligands. We aim to shed light on its plasticity analyzing the changes followed by the linker region and the relative position of the lobes using conventional molecular dynamics, accelerated MDand scaled MD (sMD). Materials & methods: Three different structures of calmodulin are compared, obtaining a total of 2.5μs of molecular dynamics, which have been analyzed using the principal component analysis and clustering methodologies. Results: sMD simulations reach conformations that conventional molecular dynamics is not able to, without compromising the stability of the protein. On the other hand, accelerated MD requires optimization of the setup parameters to be useful. Conclusion: sMD is useful to study flexible proteins, highlighting those factors that justify its promiscuity.

Steered molecular dynamic simulations of conformational lock of Cu, Zn-superoxide dismutase

The conformational lock was a bio-thermodynamic theory to explain the characteristics of interfaces in oligomeric enzymes and their effects on catalytic activity. The previous studies on superoxide dismutases (Cu, Zn-SODs) showed that the dimeric structure contributed to the high catalytic efficiency and the stability. In this study, steered molecular dynamics simulations were used firstly to study the main interactions between two subunits of Cu, Zn-SODs. The decomposition process study showed that there were not only four pairs of hydrogen bonds but also twenty-five residue pairs participating hydrophobic interactions between A and B chains of SOD, and van der Waals interactions occupied a dominant position among these residue pairs. Moreover, the residue pairs of hydrogen bonds played a major role in maintaining the protein conformation. The analysis of the energy and conformational changes in the SMD simulation showed that there were two groups (two conformational locks) between A and B chains of SOD. The first group consisted of one hydrogen-bond residues pair and seven hydrophobic interactions residues pairs with a total average energy of −30.10 KJ/mol, and the second group of three hydrogen-bond residues pair and eighteen hydrophobic interactions residues pairs formed with a total average energy of −115.23 KJ/mol.

Lessons learned about steered molecular dynamics simulations and free energy calculations.

The calculation of free energy profiles is central in understanding differential enzymatic activity, for instance, involving chemical reactions that require QM-MM tools, ligand migration, and conformational rearrangements that can be modeled using classical potentials. The use of steered molecular dynamics (sMD) together with the Jarzynski equality is a popular approach in calculating free energy profiles. Here, we first briefly review the application of the Jarzynski equality to sMD simulations, then revisit the so-called stiff-spring approximation and the consequent expectation of Gaussian work distributions and, finally, reiterate the practical utility of the second-order cumulant expansion, as it coincides with the parametric maximum-likelihood estimator in this scenario. We illustrate this procedure using simulations of CO, both in aqueous solution and in a carbon nanotube as a model system for biologically relevant nanoheterogeneous environments. We conclude the use of the second-order cumulant expansion permits the use of faster pulling velocities in sMD simulations, without introducing bias due to large dispersion in the non-equilibrium work distribution.

Effects of protein-protein interactions and ligand binding on the ion permeation in KCNQ1 potassium channel

The voltage-gated KCNQ1 potassium ion channel interacts with the type I transmembrane protein minK (KCNE1) to generate the slow delayed rectifier (IKs) current in the heart. Mutations in these transmembrane proteins have been linked with several heart-related issues, including long QT syndromes (LQTS), congenital atrial fibrillation, and short QT syndrome. Off-target interactions of several drugs with that of KCNQ1/KCNE1 ion channel complex have been known to cause fatal cardiac irregularities. Thus, KCNQ1/KCNE1 remains an important avenue for drug-design and discovery research. In this work, we present the structural and mechanistic details of potassium ion permeation through an open KCNQ1 structural model using the combined molecular dynamics and steered molecular dynamics simulations. We discuss the processes and key residues involved in the permeation of a potassium ion through the KCNQ1 ion channel, and how the ion permeation is affected by (i) the KCNQ1-KCNE1 interactions and (ii) the binding of chromanol 293B ligand and its derivatives into the complex. The results reveal that interactions between KCNQ1 with KCNE1 causes a pore constriction in the former, which in-turn forms small energetic barriers in the ion-permeation pathway. These findings correlate with the previous experimental reports that interactions of KCNE1 dramatically slows the activation of KCNQ1. Upon ligand-binding onto the complex, the energy-barriers along ion permeation path are more pronounced, as expected, therefore, requiring higher force in our steered-MD simulations. Nevertheless, pulling the ion when a weak blocker is bound to the channel does not necessitate high force in SMD. This indicates that our SMD simulations have been able to discern between strong and week blockers and reveal their influence on potassium ion permeation. The findings presented here will have some implications in understanding the potential off-target interactions of the drugs with the KCNQ1/KCNE1 channel that lead to cardiotoxic effects.

Moving Beyond the Mechanical Clamp: An Exploration into Differential Mechanical Stability of Ubiquitin Family Proteins

Mechanical properties of proteins with a β-grasp/β-sandwich fold topology have been studied extensively by both single-molecule force spectroscopy experiments and steered molecular dynamics simulations (SMD), owing to their endurance towards high rupture forces. Previous studies have shown that the backbone hydrogen bonding network holding the terminal β-strands (β-clamp) in these proteins contributes prominently to their mechanical stability. However, the role of interactions beyond the β-clamp in providing mechanical resistance is not well understood. Ubiquitin family proteins, which share a common β-grasp fold including a β-clamp, exhibit distinct mechanical stability. The low sequence homology shared by these proteins implies a differential interaction network constituting the protein fold and this could serve as a paradigm to explore the role of collective structural contributions towards mechanical stability of a protein. Here, we explore the structural origin behind differences in the mechanical properties of three members of the ubiquitin family (ubiquitin, SUMO1 and SUMO2). SMD simulations have been employed to reveal a high resolution atomistic view of protein structural changes as they unravel under the influence of a mechanical force applied along their termini. The simulations corroborate experimentally observed differences in mechanical stability of these three proteins. We monitor collective structural changes to reveal the major secondary structural elements that contribute to their mechanical stability. Differential coupling between various β-strands and the α-helix results in differences in the mechanical stability of the ubiquitin family proteins. Our study broadens the outlook of mechanical stability in proteins with a β-grasp fold to go beyond the ‘β-clamp’, hitherto considered as the sole contributor to their mechanical stability. Our study highlights the coupling of β-strands with α-helix as the next major contributor to mechanical stability of proteins of ubiquitin fold.

Identification of key residues modulating the stereoselectivity of nitrile hydratase toward rac‐mandelonitrile by semi‐rational engineering

Optically pure compounds are important in the synthesis of fine chemicals. Using directed evolution of enzymes to obtain biocatalysts that can selectively produce high-value chiral chemicals is often time-, money-, and resource-intensive; traditional semi-rational designs based on structural data and docking experiments are still limited due to the lack of accurate selection of hot-spot residues. In this study, through ligand-protein collision counts based on steered molecular dynamics simulation, we accurately identified four residues related to improving nitrile hydratase stereoselectivity toward rac-mandelonitrile (MAN). All the four selected residues had numerous collisions with rac-MAN. Five mutants significantly shifting stereoselectivity towards (S)-MAN were obtained from site-saturation mutagenesis, one of them, at position βPhe37, exhibiting efficient production of (S)-MAN with 96.8% eep , was isolated and further analyzed. The increased pulling force observed during SMD simulation was found to be in good coincidence with the formation of hydrogen bonds between (R)-MAN and residue βHis37. (R)-MAN had to break these barriers to enter the active site of nitrile hydratase and S selectivity was thus improved. The results indicated that combining steered molecular dynamics simulation with a traditional semi-rational design significantly reduced the select range of hot-spot residues for the evolution of NHase stereoselectivity, which could serve as an alternative for the modulation of enzyme stereoselectivity.

Revealing the importance of linkers in K-series oxime reactivators for tabun-inhibited AChE using quantum chemical, docking and SMD studies.

Inhibition of acetylcholinesterase (AChE) with organophosphorus compounds has a detrimental effect on human life. Oxime K203 seems to be one of the promising reactivators for tabun-inhibited AChE than (K027, K127, and K628). These reactivators differ only in the linker units between the two pyridinium rings. The conformational analyses performed with quantum chemical RHF/6-31G* level for K027, K127, K203 and K628 showed that the minimum energy conformers have different orientations of the active and peripheral pyridinium rings for these reactivator molecules. K203 with (-CH2-CH=CH-CH2-) linker unit possesses more open conformation compared to the other reactivators. Such orientation of K203 experiences favorable interaction with the surrounding residues of catalytic anionic site (CAS) and peripheral anionic site (PAS) of tabun-inhibited AChE. From the steered molecular dynamics simulations, it has been observed that the oxygen atom of the oxime group of K203 reactivator approaches nearest to the P-atom of the SUN203 (3.75Å) at lower time scales (less than ~1000ps) as compared to the other reactivators. K203 experiences less number of hydrophobic interaction with the PAS residues which is suggested to be an important factor for the efficient reactivation process. In addition, K203 crates large number of H-bonding with CAS residues SUN203, Phe295, Tyr337, Phe338 and His447. K203 barely changes its conformation during the SMD simulation process and hence the energy penalty to adopt any other conformation is minimal in this case as compared to the other reactivators. The molecular mechanics and Poisson-Boltzmann surface area binding energies obtained for the interaction of K203 inside the gorge of tabun inhibited AChE is substantially higher (-290.2kcal/mol) than the corresponding K628 reactivator (-260.4kcal/mol), which also possess unsaturated aromatic linker unit.

The effect of betaine on the foam stability: Molecular simulation

Zwitterionic betaines are widely used as foam boosters due to these can enhance the stability of foam films. In this paper, mechanistic insights of betaine to improve the stability of alkyl-polyoxyethylene carboxylate (AEC) foam are provided by molecular simulation. In the simulation, we observe the electropositive nitrogen atoms in betaine interact with the electronegative sulfur atoms, an electrostatic structure is formed at the air/water interface. Interaction energies of the mixed surfactants are calculated by the quantum chemistry methods. The calculations show betaine-AEC and betaine–betaine possess attractive interaction, and that AEC–AEC has repulsion to each other. In the other words, the repulsion between the headgroups of anionic surfactants is relaxed by betaine. Additionally, the influence of concentration of betaine on the stability of foam films is also simulated. The RDF and coordination numbers show that the electrostatic structures become denser with the increasing concentration of betaine. Therefore, entry barrier is enhanced accordingly. The SMD simulation also demonstrates the same variation tendency of entry barrier. The simulation details provide vital supplements to experiments.