An experimentally validated approach to calculate the blood-brain barrier permeability of small molecules

Drug development for the treatment of central nervous system (CNS) diseases is extremely challenging, in large part due to the difficulty in crossing the blood-brain barrier (BBB). Here we develop and experimentally validate a new in silico method to predict quantitatively the BBB permeability for small-molecule drugs. We show accurate prediction of solute permeabilities at physiological temperature using high-temperature unbiased atomic detail molecular dynamics simulations of spontaneous drug diffusion across BBB bilayers. These simulations provide atomic detail insights into the transport mechanisms, as well as converged kinetics and thermodynamics. The method is validated computationally against physiological temperature simulations for fast-diffusing compounds, as well as experimentally by direct determination of the compound permeabilities using a transwell assay as an in vitro BBB model. The overall agreement of the predicted values with both direct simulations at physiological temperatures and experimental data is excellent. This new tool has the potential to replace current semi-empirical in silico screening and in vitro permeability measurements in CNS drug discovery.


Supplementary Information
Calculating solute diffusion coefficients along the membrane normal Solute diffusivities are calculated using the method of Hummer and coworkers 1 . Each solute was placed near the free energy barrier maximum (position GB) and held fixed using a harmonic umbrella potential with a spring constant 100 kJ mol -1 nm -2 . The temperature was varied between 310K and 500K. The diffusion coefficient Dz was calculated from these simulations using an autocorrelation function: (Eq. S1) where τ was calculated from the autocorrelation function: The first 10% of the trajectory was excluded from the calculations to allow for the system to equilibrate. Variances were calculated by dividing the simulation into 3 equal parts and calculating the diffusion coefficient for each part.
Calculating lipid diffusion coefficients in the plane of the membrane When calculating the lipid diffusion coefficient, a number of methodological issues need to be addressed, including the treatment of the periodic boundary conditions as well as finitesize effects of the box. [2][3][4] Therefore, the lateral lipid diffusion coefficient was determined from the mean square displacement of phosphorus atoms with respect to unwrapped periodic boundary coordinates in the plane of the membrane (xy-plane). Since the lipid composition of the leaflet is the same on average, the lateral diffusion coefficient, was calculated by averaging the mean square displacement <r(t)> 2 observed in both bilayer leaflets. Cholesterol within the bilayer are excluded due to their susceptibility to flip-flop.
Experiments measure changes in concentration on either side of the membrane, which corresponds to measuring the net flux across the membrane (i.e. forward minus backward).
In the simulations each particle can be traced and we can obtain the backward and forward fluxes separately. If the particles do not interact, which is checked in the simulations, the concentration can be thought of to be constant on one side of the bilayer and zero on the other. The permeability is the calculated from the averaged the total forward and backward fluxes through the bilayer.
Due to pressure coupling the box volume and area of the bilayer patch, S, will vary during the simulation and need to be averaged. The flux is calculated by counting the total number of forward and backward permeation events through the bilayer, dividing by the simulation length and averaging the result: where ri (#/s) and ro are the forward and backward transport rates, respectively. The permeability is then obtained using r = P·S·C: where # are the number of transport events, t is the simulation time, and S and C are simulation averages of the membrane area and concentration (the volume changes), respectively.     Figure 1D) and can be fit to: ΔG ‡ (T) = g·T + G0, with an average R 2 of 0.98. DZ: Diffusion coefficient along the membrane normal was fit against the lipid diffusion coefficient using: ln(DZ) = ln(b) + mD · ln(DL), with parameters m and b (c.f. Figure 4). lb : Barrier width was fit against temperature using: ln(lb) = ln(a) + ml · ln(T ), with fitting parameters m and a (c.f. Figure 4). Average lateral (i.e. in-plane) lipid diffusion coefficient (D) calculated from apical lipid bilayer simulations with ethanol as a solute. There are 96 lipid molecules in the system, the composition is shown in Figure S2. The phosphorus atoms of the lipids were chosen to calculate the average meansquared displacement.       Note that the transport rate across the hydrophobic core is not necessarily equal to the total rate, as some molecules have significant barriers for partitioning out of the membrane interface (c.f. Figure S4 & Figure S6). The inset shows that the lipid diffusion coefficient DL increases rapidly with increasing temperature, while the pre-factor A of ethanol molecules decreases rapidly with increasing temperature. C. Pre-factor A for six solutes in apical hBMEC lipid bilayers calculated by an Arrhenius relation. Inverse of pre-factor, A -1 , decreases with increasing temperature; this behavior can be divided into two zones: a fast-changing region (T < 400K) and slowlychanging region (T > 400K).