Fast Estimation of the Blood–Brain Barrier Permeability by Pulling a Ligand through a Lipid Membrane

The blood–brain barrier (BBB) is a physical barrier that regulates the homeostasis of the neural microenvironment. A relative estimate of the BBB permeability, which is important for drug design, may be experimentally provided by the logBB (the blood–brain concentration ratio) and the logPS (permeability–surface-area product), while many computational methods aim to identify key properties that correlate well with these quantities. Although currently existing computational methods (e.g., quantitative structure activity relation) have made a significant contribution in screening various compounds that could potentially translocate through the BBB, they are unable to provide a physical explanation of the underlying processes and they can often be computationally demanding. Here, we use steered molecular dynamics simulation to estimate the BBB permeability of various compounds on the basis of simple lipid–membrane models by computing the nonequilibrium work, Wneq, produced by pulling the compounds through the membrane. We found that the values of Wneq correlate remarkably well with logBB and logPS for a range of compounds and different membrane types and pulling speeds, independently of the choice of force field. Moreover, our results provide insight into the role of hydrogen bonds, the energetic barriers, and the forces exerted on the ligands during their pulling. Our method is computationally easy to implement and fast. Therefore, we anticipate that it could provide a reliable prescreening tool for estimating the relative permeability of the BBB to various substances.


■ INTRODUCTION
Millions of people around the world suffer from some kind of central nervous system (CNS) disorder. 1 A well-known example of such condition is Alzheimer's disease and, as in many other cases, a possible cure may require the delivery of drugs to the brain. However, 98% of small-molecule drugs and almost all of the larger ones are unable to reach the brain because of the presence of the blood−brain barrier (BBB), which, among others, protects the brain from toxins and pathogens while allowing for the delivery of necessary nutrients and oxygen to the brain. 2,3 The properties of the BBB are mainly determined by the endothelial cells, but they are induced and maintained by complex communication with other types of cells according to the needs of the CNS. 4 In fact, the transfer of substances between the blood microcirculation system and the brain parenchyma can take place through a number of active (energy is required, e.g., ATP hydrolysis) and passive translocation mechanisms, where the latter are driven by concentration and electrostatic gradients. 5 These mechanisms and the way that they apply to different substances eventually determine the BBB permeability, which can be expressed by two commonly used experimental descriptors: the steady-state concentration of a drug in the brain over the concentration in the blood (logBB) and the permeability−surface-area product (logPS). 6 In general, the entry to the CNS depends on a range of factors characterizing a compound, such as molecular weight, lipid solubility, ability to form hydrogen bonds, amino acid composition, charge, and three-dimensional (3D) structure (conformation, flexibility, folding, and others). 7−9 To this end, the role of many of these characteristics might manifest in model membranes because the rate of passive diffusion across a membrane is proportional to the partition coefficient of the compound between the membrane and the external medium. 6,10−12 In this respect, the difficulty in traversing such membranes may constitute a measure of the barrier that a substance needs to overcome, in this way providing important information for pharmacokinetics and rational drug design, which is the motivation of the current study. Thus, model membranes provide opportunities for prescreening various compounds with the aim of reducing the cost and enhancing the speed of the BBB-permeation analysis, for example, by estimating the partitioning free energy, which has been the focus of previous computational studies. 6,10,11 In general, the blood−brain permeability can be experimentally determined by various in vivo methods, for example, by measuring the ratio between the compound concentration in brain tissues and blood samples (logBB). 1,13 These methods are very reliable, but they are laborious and relatively low throughput. Therefore, in vitro methods based on various cell assays have been employed to predict the BBB permeability. 14 Still, these methods are time-consuming and expensive rendering them practically inconvenient for screening large collections of chemicals. In contrast, various computational approaches 15,16 have already served as potential prescreening tools with the aim of reducing the cost and enhancing the speed of BBB-permeability analysis. 15 Examples of such approaches are the quantitative structure activity relation and the application of artificial intelligence algorithms (e.g., machine learning, genetic algorithms, and artificial neural networks), 3 which have elucidated the role of various parameters, such as lipophilicity, hydrogen-bonding, 17−19 and other physicochemical factors of the compounds (e.g., solubility, excess molar refraction, polarizability, hydrogen-bond acidity/basicity, and drug solvation free energy in water) 8,10,20,21 by assuming that the BBB permeability is governed by these descriptors (e.g., lipophilicity, hydrogen-bonding, polarizability, acidity/basicity, polar surface area, pK a , etc.). The results of these studies have been corroborated on the basis of the high correlation coefficient between experimental and computational predictions for representative sets of compounds. 6,10,11 Molecular simulation has also been employed to predict the permeability of the BBB by using simple membrane models and free-energy calculation techniques. For instance, Lombardo et al. have shown that the BBB partitioning in terms of logBB for a series of compounds ranging from simple solutes to histamine H 2 antagonists was correlated with the computed solvation free energy in water without taking into account specificities that pertain to the membrane structure. 10 Furthermore, this correlation provided successful predictions of blood−brain partitioning for compounds outside of the training dataset, which highlights the potential of these methods as prescreening tools. In another study, Carpenter et al. have obtained very good agreement with experimental results on the logBB and the logPS of twelve drug-like compounds using the lipid bilayer as a simple BBB model and molecular dynamics (MD) simulation of allatom models combined with free-energy methods. 6 Their predictions have correlated remarkably well with both the experimental logBB and logPS, obtaining values of the correlation coefficient close to 1. A similar approach has been adopted to express the membrane permeability of small compounds through the diffusion constant and the potential of mean force in MD simulation. 11 These studies have clearly underlined the predictive power of MD methods toward estimating energetic barriers associated with BBB permeability and identifying the physical mechanisms related to the translocation of various compounds through the BBB. Still, standard MD methods combined with free energy calculation methods are computationally expensive and often involve complex simulation protocols that require a delicate and iterative process in order to improve the accuracy of the results.
In this article, we propose a simple, reliable, and fast method to provide a relative estimate of the BBB permeability to a representative set of compounds, which also includes com-pounds considered in previous computational studies. 6 To achieve this, we pull the selected compounds through model membranes by using steered molecular dynamics (SMD) simulation and monitor the nonequilibrium work, W neq , required to cross the membrane. We show that this approach obtains robust results (e.g., independent of the membrane model or the pulling speed) that correlate remarkably well with the experimental values of logBB, and logPS, as well as the values of the topological polar surface area (TPSA). Moreover, our method provides detailed information on the forces experienced by the compounds during the translocation process and the formation of hydrogen bonds between the compounds and both the model membrane and water molecules. We anticipate that the proposed approach will have broader implications in this research area as a simple, reliable, and fast prescreening tool for predicting the relative permeability of the BBB to various compounds.

■ MATERIALS AND METHODS
Choice of Compounds. In this study, we have chosen 26 small-molecule compounds, which were selected from the literature (Tables S1 in the Supporting Information). The logBB and TPSA values are known for all compounds, while the logPS values are available in the literature only for a subset of the studied compounds. Moreover, many of them have been previously studied by standard MD simulation 6 by using the GROMOS54a7 force-field. Here, we have modeled a set of 26 small-molecule compounds not only with the GROMOS54a7 force-field, but, also, by using the CHARMM36 force-field. All the compounds have been previously used in experimental BBB permeability reports and are well studied. Hence, the 3D structures, a few common descriptors, and the logBB and logPS values of the selected compounds have been documented in the literature and obtained from the Collaborative Drug Discovery in PubChem 22,23 (Table S1 in the Supporting Information). The selected 26 compounds have logBB values that are evenly distributed between −2.51 and 1.64. Fourteen out of the 26 compounds can cross the BBB by a diffusion mechanism, while the rest cannot translocate through the BBB.
Model for the Compounds and the Lipid Bilayer. A typical system consists of a lipid bilayer surrounded by water and one of the compounds (Table S1), which is positioned in the water domain on the left of the lipid bilayer ( Figure 1). The total size of our systems was typically about 40,000 atoms. Following previous work, 6 the bilayer consists of 128 DOPC (1,2-dioleoyl- Figure 1. Typical simulation setup of our system. The system consists of a lipid bilayer (green), water (blue), and one of the compounds (orange and green). The zero of the axis in the z direction is positioned at the center of the lipid bilayer. The SMD force that pulls the dummy atom through the bilayer in the z direction is indicated with a red arrow.
Journal of Chemical Information and Modeling pubs.acs.org/jcim Article sn-glycero-3-phosphocholine) lipids 24,25 modeled with the Berger force-field. 24−27 The so-called "Berger lipids" can be combined with the GROMOS representation of the compounds. 27 Hence, the GROMOS54a7 force-field was used to model the 26 compounds. Moreover, the CHARMM36 28−30 force-field was employed for both the bilayer by using the CHARMM-GUI server 31−33 and the 26 compounds (Table S1). The parameters of the compounds were obtained from the ATB server 34−37 for the GROMOS54a7 force-field and the SwissParam server 38 in the case of the CHARMM36 forcefield. In general, the ATB server combines a knowledge-based approach with quantum mechanical calculations. First, the ligand was optimized using the HF/STO-3G theory and then reoptimized at a more advanced level B3LYP/6-31G*. 39−41 To estimate the point charges of compound atoms, the electrostatic potential was fitted using the Kollman−Singh procedure. 42 The bond and angle force constants were extracted from the Hessian matrix. All ligands were modeled in their neutral form.
We have used a lipid bilayer with POPC lipids (2-oleoyl-1pamlitoyl-sn-glyecro-3-phosphocholine) in addition to DOPC in order to assess the effect of the membrane composition on our conclusions regarding the translocation ability of the compounds. The model for POPC is also based on the CHARMM36 force-field and consists of 128 POPC molecules. As in the case of the DOPC, the membrane topology for the POPC membrane model was created by using the CHARMM-GUI server, 31,33,50 while the topology and the parameters for small compounds were obtained by the SwissParam server, 38 which is compatible with the CHARMM36 all-atom force-field. Also, a DOPC membrane model with cholesterol was studied with a DOPC/ cholesterol ratio of 2:1 and using the CHARMM force-field. Finally, the simple point charge water model was used in the case of the GROMOS54a7 force-field, while the TIP3P water model was employed in the case of the CHARMM36 force-field. 43 Steered Molecular Dynamics. SMD offers the possibility of obtaining results much quicker than it might be possible with standard MD simulation. 44 In the case of the SMD simulation, a spring is attached from one side to a dummy atom and from the other side to the compound. Then, the dummy atom is pulled from its initial position toward the bilayer in the z direction with a constant velocity v, covering a distance Δz c = vt ( Figure 1) with t indicating time. Hence, the elastic force experienced by the compound at a certain time during pulling is F = k(Δz − vt), where Δz is the displacement of the compound's atom connected with the spring in the direction of pulling. As in the case of atomic force microscopy experiments, 45 we chose k = 600 kJ/(mol·nm 2 ) and different pulling velocities (i.e., v = 0.5, 1.5, and 5.0 nm/ns). These values have been also used in our previous studies and are well-tested. 46−48 Moreover, the spring constant is suitable to carry out the pulling experiments, as can be verified by the average position−time profiles for the compounds (see Figure S1 in Supporting Information). To prevent the membrane from drifting because of the external force, the phosphate atoms on both lipid layers were restrained by a harmonic potential with a spring constant of 1000 kJ/(mol· nm 2 ) in the z direction.
The main output of the simulation is the force, F, as a function of time and the position, z, of the compound, from which we can readily calculate the work, W neq , required to pull the compound through the bilayer over a distance L z = ∫ dz, where L z is the width of the bilayer in the z direction. W neq can be mathematically expressed as W neq = ∫ F dz. The average force−time and force−position profiles obtained from five individual trajectories for the 26 studied compounds are shown in Figure S2 in the Supporting Information in the case of the GROMOS54a7 force-field and the DOPC membrane model.
Once the system was set up as shown in Figure 1, the steepest descent and conjugate gradient methods were used for energy minimization, as is usually done. The energy-minimized system was then utilized as an initial configuration for the SMD simulation. The simulation was carried out at T = 323 K, by employing the Nose−Hoover thermostat with an effective relaxation time parameter set to τ T = 0.5 ps. The pressure was maintained at 1 bar using the semi-isotropic Parrinello−Rahman barostat with relaxation time τ P = 1 ps, and a compressibility of 4.5 × 10 −5 bar −1 . The LINCS algorithm 49,50 was used to constrain bond lengths, thus allowing for a 2 fs time-step during simulation. A cutoff of 1.4 nm was chosen for nonbonded interactions, and the neighbor list was updated every 10 ps, while long-range electrostatic interactions were calculated using the particle−mesh Ewald method. Here, the GROMACS 5.1 package was used to carry out all simulations. Simulations Journal of Chemical Information and Modeling pubs.acs.org/jcim Article typically last from 1.8 to 18.0 ns depending on the pulling velocity. In our study, the considered range of pulling velocities span an order of magnitude, namely, v = 0.5, 1.5, and 5.0 nm/ns. Average properties were calculated from an ensemble of five trajectories with different initial conditions by using a random initial velocity generation and changing the orientation of the pulled compound ( Figure S3). Examples of individual trajectories are shown in Figures S4 and S5 in the Supporting Information as well as Tables S2 and S3 for different membrane models and force fields in the case of mannitol, which has the smallest experimental logBB and logPS values and the highest TPSA value. In fact, the logBB and the TPSA values for the 26 compounds of our study appear to be anticorrelated with R 2 = 0.53 (R is the correlation coefficient), as shown in Figure S6.
Hence, the available values in the literature of logBB and TPSA for our compounds suggest that a larger polar surface area of the compounds is associated with a lower permeability of the BBB.
Moreover, values of TPSA below 90 Å 2 suggest an increased potential for BBB penetration, while larger values are usually associated with lower probability of crossing the BBB. 51

■ RESULTS AND DISCUSSION
Correlation of the maximum force, F max , and the nonequilibrium work, W neq , with the experimental logBB and logPS values, as well as the TPSA values. Figure 2 shows typical snapshots that illustrate the translocation of a compound through the model membrane, in this case mannitol (see also movie in the Supporting Information),  (Table S1) in the case of the CHARMM36 force-field for the DOPC membrane model and pulling velocity v = 0.5 nm/ns (Table S4). The correlation coefficient R 2 is indicated at the top right corner of each graph. Figure 4. The correlation between F max and the nonequilibrium work, W neq , with logBB, logPS, and TPSA for the compounds (Tables S1 and S6) in the case of the GROMOS54a7 force-field for the DOPC membrane model and pulling velocity v = 0.5 nm/ns (Table S5). The correlation coefficient R 2 is indicated at the top of each graph. Note that logPS is only available for 8 compounds.
Journal of Chemical Information and Modeling pubs.acs.org/jcim Article where the pulling velocity v = 0.5 nm/ns and the GROMOS54a7 force-field for the model DOPC was used in this case. Obviously, the membrane remains stable during the simulation. Our analysis for each case of force field, membrane model, and pulling velocity, yields the maximum force, F max , exerted on the compound, and the nonequilibrium work, W neq , required to pull the compound through the membrane can be readily calculated from the force−position profiles (e.g., Figure S2 in the Supporting Information) as discussed in our previous section. The values of F max and W neq are reported in Table S4 of the Supporting Information and plotted as a function of logBB in Figure 3 for the 26 compounds and the DOPC membrane  (Table S6) in the case of the CHARMM36 force-field, velocity v = 5 nm/ns, and the POPC membrane model. Journal of Chemical Information and Modeling pubs.acs.org/jcim Article model simulated with the CHARMM36 forcefield and pulling velocity v = 0.5 nm/ns. Our results indicate that the correlation between F max and the logBB values of the compounds is high, namely, R 2 = 0.85. In particular, the maximum force increases as the logBB values decrease. This result agrees well with the expectation that a larger force may be required to translocate a compound through the membrane when logBB is lower (e.g., mannitol). In the case of W neq , the correlation of the simulation results with the experimental logBB values of the compounds is similarly high, namely, the correlation coefficient is R 2 = 0.86. Hence, our results indicate a remarkable correlation between the experimentally measured logBB parameters and the reported maximum force values, F max , and the nonequilibrium work, W neq , despite the simplicity of our approach. Both quantities are larger when logBB is smaller. Our analysis in the case of 26 compounds modeled with the GROMOS54a7 force-field and using the same pulling velocity (v = 0.5 nm/ns) and the DOPC membrane model as in the case of the CHARMM36 model is shown in Figure 4 and the corresponding values of the maximum force, F max , and the nonequilibrium work, W neq , as a function of the experimental logBB and logPS values as well as the TPSA are reported in Table S5 in the Supporting Information. Overall, our results suggest a high correlation between the simulation and the experimental results. In all cases, the correlation coefficient, R, obtains large values for both F max and W neq . In particular, R 2 obtains values larger than 0.80 for the correlation of F max and W neq with logBB. The correlation with the logPS is of the same magnitude but slightly smaller. Finally, the correlation of F max and W neq with the TPSA of the compounds is also significant. While F max shows a correlation as large as in the case of logBB and logPS, W neq manifests a smaller correlation with the TPSA values, that is, R 2 = 0.58. This is expected given the correlation between the logBB and TPSA values ( Figure S6 in Supporting Information). Our study also suggests that the use of the Berger lipid force-field for the membrane model, which may be insufficient in some cases (e.g., when a cutoff of 1.0 nm is used), 52,53 yields consistent results with those obtained by using the CHARMM36 force-field (cf. Figure 3). Hence, our conclusions, which are validated for a large set of different compounds, are not affected by the choice of the force field for both the compounds and the membrane. Moreover, our results are in very good agreement with the results obtained in the study of Carpenter et al., 6 where twelve compounds were investigated by means of the potential of mean force between the compounds and the BBB. In summary, our results indicate a remarkable correlation between the simulation and the experimental results  Journal of Chemical Information and Modeling pubs.acs.org/jcim Article irrespective of the force field used to model the membrane and the compounds, which may justify the use of our approach as a method for prescreening compounds with the aim of estimating eventually the relative BBB permeability to different compounds. Effect of Pulling Velocity and Membrane Type. In order to assess the influence of the pulling velocity, we carried out SMD simulations with a range of pulling velocities, that is, v = 0.5, 1.5, and 5 nm/ns, and different bilayer compositions (DOPC or POPC) for the 26 compounds of this study by using different force fields. Our results in the case of the GROMOS54a7 force-field for the DOPC membrane model (cf. Figure 4) indicate that the correlation remains high when the pulling velocity increases ( Figure S7 and Table S5), namely v = 1.5 nm/ns. Hence, the values of R 2 for the correlation of F max and W neq with logBB are not affected significantly by changing the pulling velocity. Similarly, the correlation of F max and W neq with the logPS and TPSA is of the same order. Further increase of the velocity (i.e., v = 5.0 nm/ns) suggests that the correlation of F max and W neq with logBB, logPS, and TPSA remains again of the same order ( Figure S8 and Table S5). Hence, changes in the pulling velocity do not affect our conclusions. Furthermore, we have verified our conclusions for both the GROMOS54a7 and the CHARMM36 force-fields, as well as different membrane models (DOPC and POPC). Then, Figure 5 shows the correlation of F max and W neq with the experimental logBB and logPS values in the case of the CHARMM36 force-field and pulling velocity v = 5 nm/ns for the POPC membrane model. The values are also reported in Table S6. We find that the correlation of F max and W neq with logBB, logPS, and TPSA is similar, irrespective of the membrane type (DOPC or POPC). In particular, the maximum force and the nonequilibrium work remains highly correlated with the experimental logBB values, namely, R 2 = 0.72 and R 2 = 0.80, respectively. The correlation with the logPS is also high in the case of F max , but it is much lower in the case of W neq versus logPS. The results for F max and W neq with the TPSA values show a much weaker correlation in comparison with the experimental logBB and logPS values. In the case of Figure 5, R 2 = 0.47 and 0.43 for the correlation of F max and W neq with TPSA, respectively. In summary, we conclude that the correlation between F max and W neq with logBB, logPS, and TPSA are robust because it is not affected by the employed forcefield for the compound and the membrane, the pulling velocity and the membrane model (POPC or DOPC). In all cases, our results suggest a significant correlation between F max and W neq with the experimental logBB and logPS values, as well as TPSA values, reported in the literature.
Hydrogen Bonds. We have analyzed the hydrogen bond formation between the compounds and water ( Figure S9 in the Supporting Information), as well as between the compounds and the lipids ( Figure S10 in the Supporting Information), as a function of the compound position. A hydrogen bond is formed when the distance between the donor D and acceptor A is ≤ 3.5 Å and the D−H−A angle is ≥135°(H denotes the hydrogen atom). The number of hydrogen bonds was estimated as the average number of hydrogen bonds from the ensemble of our trajectories separately for each compound. Analysis of the hydrogen-bond formation between compounds and water molecules ( Figure S9) suggests that hydrogen bonds are present in all cases, but compounds with lower logBB (smaller  (Table S1) in the case of the CHARMM36 force-field for the DOPC/cholesterol (2:1) membrane model and pulling velocity v = 5 nm/ns (Table S4). The correlation coefficient R 2 is indicated at the top right corner of each graph.
Journal of Chemical Information and Modeling pubs.acs.org/jcim Article permeability) overall have higher tendency of forming hydrogen bonds with water. Moreover, the profiles are rather symmetric with water molecules being also present in the middle of the lipid membrane during simulations. The latter effect may be a consequence of water molecules being dragged by the compounds during the pulling process, which mostly takes place for compounds with low logBB (high tendency of forming hydrogen bonds). In the case of hydrogen-bond formation between the compound and the lipids, we can observe that the presence of hydrogen bonds is related to the barrier that a compound needs to overcome in order to translocate through the membrane. In particular, the hydrogen bonds formed between mannitol and the lipid molecules is the highest among all compounds reaching the value of 3.8 ( Figure S10). Moreover, because of the apparent symmetric composition of the two membrane layers in lipids, the profiles seem to be symmetric with a first high peak in the number of hydrogen bonds appearing in the left layer ( Figure S10) and a subsequent lower peak. Finally, compounds that translocate easily through the membrane (high logBB values) do not form as many hydrogen bonds.
We further attempt to measure the correlation between the number of hydrogen bonds and the parameters logBB, logPS, and TPSA as we did previously in the case of F max and W neq . Although in the latter cases we had found a remarkable correlation, the correlation between the number of hydrogen bonds and the experimental logBB and logPS seems weaker (e.g., see Figure 6, where the case of GROMOS54a7 force-field for the DOPC membrane model and pulling velocity v = 0.5 nm/ns is illustrated. Similar results we have obtained for all cases considered in our study). This is clearly manifested by the values of the correlation coefficient (R 2 = 0.24 and 0.42, for the number of hydrogen bonds vs logBB and logPS, respectively; Figure 6), which however clearly indicate that the number of hydrogen bonds decreases as the values of logBB and logPS increase. On the contrary, the number of hydrogen bonds increases with the increase of the TPSA with the correlation coefficient being R 2 = 0.48, which is consistent with the correlation observed for the F max and W neq . Our results show that compounds with low logBB and logPS and high TPSA are able to form a larger number of hydrogen bonds with the lipids, which may suggest that the low BBB permeability be related to the ability of hydrogen-bond formation between the compounds and the lipids or the solvent. However, our results also suggest that this correlation be rather weaker than the one observed between F max and W neq , and logBB and logPS, which suggests that other factors might contribute to this correlation.
Electrostatic and Van der Waals Interactions. SMD simulation allows for the analysis of the interactions between the compounds and the lipid bilayer into individual contributions, for example, van der Waals (vdW) and electrostatic interactions. Such interaction profiles with respect to the compound's position are provided in the Supporting Information ( Figure  S11). A careful look at these profiles suggests that the vdW interactions (short range) are more pronounced in the central region of the bilayer with a width of about 4 nm. In contrast, the electrostatic interactions (long range) between the compound and the lipids are less pronounced in the very central region between the two layers and within a region with width of about 2 nm in the z direction. By considering the magnitude of both the vdW and the electrostatic interactions for every position within the bilayer, we observe that these interactions are overall most pronounced at a distance ±2 nm from the center of the lipid bilayer. However, the electrostatic energy profiles appear to be slightly asymmetric with the energy obtaining lower values in the left part of the lipid bilayer. In contrast, the vdW interactions indicate some asymmetry in the right part of the bilayer, which may indicate the effect of steric interactions between the compounds and the lipids as the compound enters the bilayer.
In Figure 7, we plot the average values of the vdW and electrostatic energies as a function of logBB, logPS, and TPSA and analyze their correlation on the basis of the correlation coefficient R, as we have done previously. While Figure 7 shows results for the GROMOS54a7 force-field and the DOPC membrane model for pulling velocity v = 0.5 nm/ns, our conclusions are similar for all cases of our study. In the case of the vdW interactions, we observe that the magnitude of these interactions barely correlates with the logBB and TPSA, while a more pronounced correlation is observed with the logPS values, that is, R 2 = 0.32. In particular, the larger the logPS, the greater the effect of the van der Waals interactions during the pulling of the compound. A higher degree of correlation is observed in the case of electrostatic interactions with respect to the experimental logBB and logPS as well as the TPSA. In this case, the correlation coefficient, R, obtains positive values, while R 2 = 0.49 (logBB) and R 2 = 0.70 (logPS). We find that larger values of logBB and logPS (higher permeability) correspond to smaller (absolute) values of electrostatic energies. Hence, a smaller effect of electrostatic forces is linked to a larger permeability of the BBB. On the contrary, as the TPSA increases electrostatic forces become more dominant (their absolute value increases) with R 2 = 0.43 (a small opposite effect is observed for the vdW forces). In summary, our results indicate a larger role of electrostatic interactions during the compound translocation in contrast to the van der Waals interaction. In the case of logBB and TPSA, the magnitude of correlation is of the same order as in the case of the hydrogen bonds. This may suggest that electrostatic interactions and hydrogen bonds play a more important role in the translocation of compounds than van der Waals interactions.
Effect of Cholesterol. To assess further the effect of the membrane composition, we have altered the membrane composition by including cholesterol (Figure 8) with a DOPC/cholesterol ratio of 2:1. Here, results of the maximum force and the nonequilibrium work in the case of pulling velocity v = 5 nm/ns and the CHARMM36 force-field as a function of the logBB and logPS experimental values are shown in Figure 9. The estimated R 2 values for the correlation are in very good agreement with the membrane cases without the cholesterol. Indeed, the addition of cholesterol does not affect the correlation. The results of Figure 9 are documented in Table  S4 for each compound. While the correlation of the maximum force, F max , and the nonequilibrium work, W neq , with the logBB and logPS experimental values is not affected by the membrane model, the absolute values of F max (as a result also those of W neq ) have increased by a factor of two (an example of a force−time profile is shown in Figure S12 in the Supporting Information). Hence, we conclude that the presence of cholesterol induces a larger resistance to the compound. In other words, the permittivity of the membrane is significantly reduced in a proportional way without affecting the correlation between the experimental and the simulation predictions. structure, and so on. Given that most of the drugs are unable to cross the BBB, a large number of substances may need to be tested in order to identify the ones that could translocate through the barrier. However, many of the current experimental approaches are low throughput and impractical for screening large sets of substances. Moreover, existing computational approaches can be computationally expensive or may neglect the physical mechanisms of the underlying translocation processes. Here, we have shown that the SMD simulation based on allatom models may be a fast, reliable, and computationally inexpensive approach to determine the relative permeability of the BBB to various substances. To achieve this, simple model membranes were used. Our approach based on SMD was able to provide good agreement with the experimental results in terms of the relative permeability of the compounds. In particular, we have clearly shown that compounds with lower logBB are highly correlated with larger values of nonequilibrium work required to pull a compound through the bilayer, with correlation coefficient values close to unity. Moreover, the maximum force has also shown a remarkable correlation with the experimental logBB values. The correlation between the simulation and experimental data is consistent and parameters such as the forcefield, the number of independent trajectories, and orientation of compounds ( Figure S3 in the Supporting Information), the number of compounds, the membrane model and composition, and the pulling velocity, do not affect our conclusions. We have also monitored the formation of hydrogen bonds between the compounds and both the lipids and the water molecules, where a weaker correlation between simulation and experiment was observed. Overall, the presence of a large number of hydrogen bonds was associated with experimentally low permeability of the BBB. Analysis of the van der Waals and electrostatic interactions between the compounds and the lipid bilayer have indicated a higher correlation between logBB/logPS and electrostatic interactions rather than between logBB/logPS and van der Waals interactions. Because all ligands of our study were in neutral form, the effect of the charge of the compounds 54 on the membrane permeability can be further considered in future studies. Overall, we anticipate that our method opens new opportunities for estimating the relative permeability of the BBB to various compounds and future studies may involve more complex membranes and coarse-grained models, including models for the BBB.
List of all simulated compounds (CID), their logBB, logPS, TPSA, 3D structure, and selected descriptors, including references related to the compounds; values of F, F max , and W eq for mannitol for five individual trajectories; independent force profiles for mannitol from different trajectories and their average position− time profile; force−time profiles for compounds of this study; W neq and F max versus logBB and logPS for different velocities in the case of 26 compounds simulated with the GROMOS54a7 force-field; number of hydrogen bonds versus time and position and electrostatic interaction and vdW interaction energy versus position for the studied compounds; and results with different membranes including cholesterol and an example of force−time profile (PDF) Translocation of mannitol through the DOPC membrane (cf. Figure 2)