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, pKa, 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 H2 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 all-atom 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 compounds 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, Wneq, 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 ). 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-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 force-field. 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.

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.

We have used a lipid bilayer with POPC lipids (2-oleoyl-1-pamlitoyl-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 Δzc = vt (Figure ) 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·nm2) 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·nm2) 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, Wneq, required to pull the compound through the bilayer over a distance L = ∫dz, where L is the width of the bilayer in the z direction. Wneq can be mathematically expressed as Wneq = ∫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 , 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 Nosé–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 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 R2 = 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, Fmax, and the nonequilibrium work, Wneq, with the experimental logBB and logPS values, as well as the TPSA values.

Figure 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), 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, Fmax, exerted on the compound, and the nonequilibrium work, Wneq, 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 Fmaxand Wneq are reported in Table S4 of the Supporting Information and plotted as a function of logBB in Figure for the 26 compounds and the DOPC membrane model simulated with the CHARMM36 forcefield and pulling velocity v = 0.5 nm/ns. Our results indicate that the correlation between Fmax and the logBB values of the compounds is high, namely, R2 = 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 Wneq, the correlation of the simulation results with the experimental logBB values of the compounds is similarly high, namely, the correlation coefficient is R2 = 0.86. Hence, our results indicate a remarkable correlation between the experimentally measured logBB parameters and the reported maximum force values, Fmax, and the nonequilibrium work, Wneq, despite the simplicity of our approach. Both quantities are larger when logBB is smaller.

Figure 2

Typical snapshots obtained from simulation at different times as indicated during the pulling of mannitol with velocity v = 0.5 nm/ns in the case of the GROMOS54a7 force-field and DOPC membrane model. Phosphate atoms are restrained in order to prevent the displacement of the whole membrane as the compound penetrates the bilayer.

Figure 3

Correlation of Fmax (left panel) and nonequilibrium work (right panel), Wneq, with logBB for the 26 compounds considered in this study (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 R2 is indicated at the top right corner of each graph.

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 and the corresponding values of the maximum force, Fmax, and the nonequilibrium work, Wneq, 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 Fmax and Wneq. In particular, R2 obtains values larger than 0.80 for the correlation of Fmax and Wneq with logBB. The correlation with the logPS is of the same magnitude but slightly smaller. Finally, the correlation of Fmax and Wneq with the TPSA of the compounds is also significant. While Fmax shows a correlation as large as in the case of logBB and logPS, Wneq manifests a smaller correlation with the TPSA values, that is, R2 = 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 ). 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 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.

Figure 4

The correlation between Fmax and the nonequilibrium work, Wneq, 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 R2 is indicated at the top of each graph. Note that logPS is only available for 8 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 ) 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 R2 for the correlation of Fmax and Wneq with logBB are not affected significantly by changing the pulling velocity. Similarly, the correlation of Fmax and Wneq 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 Fmax and Wneq 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 shows the correlation of Fmax and Wneq 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 Fmax and Wneq 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, R2 = 0.72 and R2 = 0.80, respectively. The correlation with the logPS is also high in the case of Fmax, but it is much lower in the case of Wneqversus logPS. The results for Fmax and Wneq with the TPSA values show a much weaker correlation in comparison with the experimental logBB and logPS values. In the case of Figure , R2 = 0.47 and 0.43 for the correlation of Fmax and Wneq with TPSA, respectively. In summary, we conclude that the correlation between Fmax and Wneq 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 Fmax and Wneq with the experimental logBB and logPS values, as well as TPSA values, reported in the literature.

Figure 5

Correlation between Fmax and Wneq with logBB (top panel) and logPS (bottom panel) for the compounds considered in this study (Table S6) in the case of the CHARMM36 force-field, velocity v = 5 nm/ns, and the POPC membrane model.

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 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 Fmax and Wneq. 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 , 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 (R2 = 0.24 and 0.42, for the number of hydrogen bonds vs logBB and logPS, respectively; Figure ), 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 R2 = 0.48, which is consistent with the correlation observed for the Fmax and Wneq. 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 Fmax and Wneq, and logBB and logPS, which suggests that other factors might contribute to this correlation.

Figure 6

Number of hydrogen bonds at the position of the maximum force experienced by the compound in the case of the GROMOS54a7 force-field with the DOPC model and pulling velocity v = 0.5 nm/ns as a function of logBB, TPSA (top panel), and logPS (bottom panel). The square of the correlation coefficient, R, for each case is indicated on the top left corner of each graph.

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 , 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 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, R2= 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 R2 = 0.49 (logBB) and R2 = 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 R2 = 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.

Figure 7

Correlation of the average vdW or electrostatic interaction energies with logBB, TPSA, and logPS at the position of the maximum force experienced by the compound in the case of the GROMOS54a7 force-field, the DOPC membrane model and pulling velocity v = 0.5 nm/ns. The correlation coefficient R for each case is indicated on the top left corner of each graph.

Effect of Cholesterol

To assess further the effect of the membrane composition, we have altered the membrane composition by including cholesterol (Figure ) 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 . The estimated R2 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 are documented in Table S4 for each compound. While the correlation of the maximum force, Fmax, and the nonequilibrium work, Wneq, with the logBB and logPS experimental values is not affected by the membrane model, the absolute values of Fmax (as a result also those of Wneq) 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.

Figure 8

Snapshot of an equilibrium membrane of DOPC/cholesterol with a ratio of 2:1. Cholesterol molecules are highlighted with vdW spheres. Pymol software was used for the visualization of the molecules.

Figure 9

Correlation of Fmax (left panel) and nonequilibrium work (right panel), Wneq, with logBB, TPSA, and logPS for compounds considered in this study (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 R2 is indicated at the top right corner of each graph.

Conclusions

The ability of drugs to cross the BBB depends on their physicochemical characteristics, such as lipophilicity, 3D 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 all-atom 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.