Influence of local ordering in the permeation of Temozolomide through the brain plasmatic membrane

Influence of local ordering in the permeation of Temozolomide through the brain plasmatic membrane Yanhong Ge a , Huixia Lu b,* , Jordi Martí a,* a Department of Physics, Polytechnic University of Catalonia-Barcelona Tech, B4-B5 Northern Campus UPC, Barcelona 08034, Catalonia, Spain b Institut de Ciencia de Materials de Barcelona (ICMAB-CSIC), Campus de la UAB, Bellaterra 08193, Catalonia, Spain ARTICLE INFO Keywords: Temozolomide Brain plasma membrane Free energy barrier Molecular dynamics Adaptive biasing force ABSTRACT Temozolomide, a small-molecule drug, is primarily used to treat glioblastoma, a tumor that attacks both the spinal cord and brain. Understanding how Temozolomide interacts with different lipids within the brain cell membrane at the atomic level can help elucidate its ability to permeate through cell membranes. In this study, we constructed a simplified brain plasma membrane model to explore the microscopic structure and dynamics of Temozolomide using all-atom microsecond-scale molecular dynamics simulations. Temozolomide is typically found in the solvent-aqueous fluid surrounding the brain membrane, but it can access the membrane interface regularly and eventually bind to lipids of the choline and cerebroside classes. To investigate the free energy barriers of Temozolomide related to its crossing of brain-like plasma membranes, we employed adaptive biasing force methods. These simulations revealed that the free energy barriers ranged between 28 and 50 kcal/mol at temperatures between 310 K and 323 K. Our findings suggest that Temozolomide cannot cross the membrane by pure diffusion at normal human body temperature, but that rising the temperature significantly increases the probability of barrier crossing. This is primarily due to the crucial role played by cholesterol and lipids of the cerebroside class. These results can be used to optimise the molecular design of Temozolomide and develop new analogs with improved pharmacokinetic properties. 1. Introduction Temozolomide (TMZ) is the most well-known and effective drug commonly used to treat patients with glioblastoma, a highly malignant and difficult-to-treat brain tumor [1,2], as well as anaplastic astrocytoma [3]. TMZ is a small-molecule imidazotetrazine that acts as an alkylating agent with demonstrated efficacy against glioblastoma [4]. However, several barriers —most notably the reticuloendothelial system and the blood-brain barrier (BBB) [5]— limit the ability of many drugs to reach the brain and effectively target tumors. Overcoming these obstacles is essential for successful treatment [6–8]. Under normal conditions, TMZ, which has a short plasma half-life (∼1.8 h), is completely absorbed in the gastrointestinal tract following oral administration [9]. Subsequently, the prodrug rapidly hydrolyzes at physiological pH into its active but unstable metabolite, 5-(3-methyltriazen-1-yl)imidazole-4carboxamide (MTIC, half-life ∼2 min) [9–11], which then degrades into 5-aminoimidazole-4-carboxamide and the methyldiazonium ion [12]. In the final step, the methyldiazonium cation preferentially methylates DNA at guanine and adenine residues [13,14]. It is important to note that the activation of TMZ occurs within a narrow pH range close to physiological pH [15]. The MTIC compound formed in the plasma is unable to cross the BBB and must instead be generated locally within the brain. However, TMZ itself is generally considered to have good BBB permeability [16,17]. Despite its effectiveness, prolonged TMZ therapy can result in severe and unpredictable myelosuppression, which has limited its further clinical development [18]. To enhance TMZ concentrations in the brain, many researchers have explored nanoscale drug delivery systems as a promising strategy [11]. For example, Bouzinab et al. [19] developed a nanodelivery platform based on a biocompatible protein nanocage that improved TMZ delivery to cancer cells, demonstrating increased efficacy of the imidazotetrazine class for treating glioblastoma multiforme and potentially other malignancies. Meanwhile, other research groups are focusing on the development of novel compounds, including TMZ derivatives [20–22]. Recently, Yin et al. [23] reported a new TMZ-derived compound, 5-aminoimidazole-4-carboxamide (AICA). The efficacy of TMZ is also limited by its restricted ability to cross brain and tumor cell membranes to reach its target site [24]. Numerous * Corresponding authors. E-mail addresses: [email protected] (H. Lu), [email protected] (J. Martí). Contents lists available at ScienceDirect Biophysical Chemistry journal homepage: www.elsevier.com/locate/biophyschem https://doi.org/10.1016/j.bpc.2025.107457 Received 23 January 2025; Received in revised form 8 May 2025; Accepted 10 May 2025 Biophysical Chemistry 324 (2025) 107457 Available online 26 May 2025 0301-4622/© 2025 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/bync-nd/4.0/ ). studies have investigated how factors such as age and regional brain distribution influence the lipid composition of neurons and brain cells [25–28]. The human brain is composed of specialized cells with membranes structured in multiple layers, each containing distinct lipid components [29]. TMZ is primarily used in the treatment of glioblastoma, which most commonly develops in the supratentorial region of the brain [30]. Fluorescence and spectrophotometry experiments have shown that TMZ has a high affinity for brain membranes, which may hinder its penetration due to strong interactions with lipid head groups located at the blood-brain barrier interface [9,24]. However, to the best of our knowledge, no studies have yet examined the structural interactions between TMZ and biological brain plasma membranes. In this work, we consider a realistically complex lipid model of the human neuronal/brain plasma membrane (BPM), which accurately represents the major components of the actual system. A schematic representation of the membranes forming the outer layers of brain cells is shown in Fig. 1. It is important to note that, in general, the inner and outer leaflets of the BPM differ in their lipid composition. In a pioneering study on computational modeling of the brain plasma membrane, Ing´ olfsson et al. [31] developed a coarse-grained (CG) model for the BPM, incorporating tens of lipids from various classes, based on lipidomic data from neurons and brain tissue. In the present work, we have drawn upon the lipid types identified in Ing´ olfsson et al.’s study [31], selecting the most abundant and representative lipid from each class. This approach has enabled us to perform all-atom level simulations while avoiding the need to model an impractically large number of atoms. In this work, we have focused our efforts on the study of the free energy barriers involved in the translocation of TMZ through the external layer of the BPM. In summary, the present study is first devoted to the analysis of the preferential sites of TMZ at the aqueous solution (representing blood plasma), and at the interface of the brain cell membrane. Second, the main result of this paper is the calculation of the free energy barriers that TMZ should pass through the BPM as modelled in a realistic way. The work is organised as follows: methods and computational details are described with all details in Section “Method”; the main results and the corresponding discussion are described in Section “Results and Discussion” and the main conclusions are reported in Section “Conclusions”. Finally, some key complementary properties are reported in the “Supplementary Data” (SD). 2. Methods We obtained the structure of TMZ (PubChem CID: 5394) from the National Center for Biotechnology Information PubChem database (http ://www.ncbi.nlm.nih.gov/pccompound) [9]. A realistic model of the brain cell membrane has been generated with the well-known CHARMM-GUI web-based tool [32,33]. Sketches of the backbone structures of TMZ and the main components of the brain membrane are represented in Fig. 2. The model membrane used in this work is a simplified but reliable representation of the mammalian brain plasma membrane, comprising two leaflets (outer and inner layers) with different compositions. A pioneering computational modeling based on CG simulations [31] found significant differences between different cell membranes, highlighting the composition of those participating in the brain membrane. Generally, cell membranes are formed by a wide diversity of lipids. In the case of the brain, its external membrane is formed by a mixture of diverse lipids, sterols, and other organic substances, such as: phosphatidylcholine (PC) lipids, phosphatidylethanolamine (PE) lipids, sphingomyelin (SM) lipids, phosphatidylserine (PS) lipids, glycolipids (GM), cerebrosides, phosphatidylinositol (PI) lipids, phosphatic acids (PA), phosphatidylinositol phosphates (PIP), ceramides (CER), lysophosphatidylcholine (LPC) lipids, lysophosphatidylethanolamine (LPE) lipids, and the sterol diacylglycerol (DAG). Besides, cholesterol (CHOL), is central in all membranes [34], constituting almost half of all lipids in the brain membrane. We considered the model proposed by Ing´ olfsson et al. [31] to analyze the interactions of lipids with TMZ. To simplify the study, we selected the most abundant lipid to represent one lipid type, neglecting those with minimal percentages. Subsequently, we constructed a 400lipid membrane based on the percentages indicated in Table 1. It is important to note that we employed all-atom molecular dynamics (MD) simulations, providing accuracy at the atomic-level. The NAMD2 software package [35] with the reparameterized CHARMM-36 m force field [36–39] was used in all MD simulations at the temperatures of 310, 316.5 and 323 K. TMZ and lipids are fully solvated by 53,556 TIP3P water molecules and sodium chloride at 0.15 Fig. 1. TMZ crossing brain cell membranes. Schematic diagram of the process of TMZ crossing the BPM. The brain bilayer cell membrane contains a wide variety of lipids, with the types and proportions of lipids in the upper and lower layers showing some significant differences. Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 2 M concentration, yielding a system of 95,235 atoms. The initial system size is 101Å ×101Å ×100Å, and periodic boundary conditions are all applied in X, Y, and Z directions. The XY plane is taken along the surface of the BPM, whereas Z is the direction normal to the instantaneous XY plane. The Van der Waals interactions included a smoothing function starting at 10Å and its cutoff was 12Å, with a pairlist distance of 16Å. Long-range electrostatic forces were computed using particle mesh Ewald [40], with a grid space of 1Å. Periodic boundary conditions in all spatial directions are considered. A Langevin thermostat [41] was used with a damping coefficient of 1 ps−1 to control the temperature, and the pressure was set at 1 atm and regulated by a Nose-Hoover Langevin piston [42] with Langevin dynamics [43] at an oscillation period of 50 fs. After 100 ns equilibration periods at the NVT ensemble, 0.5 μ s trajectories with a timestep of 2 fs were generated at the NPT ensemble for each system. All bonds involving hydrogens were set to fixed lengths, allowing fluctuations of bond distances and angles for the remaining atoms. 3. Results and discussion 3.1. Physical characteristics of BPM and TMZ at its interface To describe the phase states of the BPM simulated in the present work, we have computed the so-called deuterium order parameter SCD, as it was defined in references [44–45], which can help us to efficiently characterise the ordering inside the hydrated lipid bilayer. SCD is defined for each CH2 group of the two tails of each lipid class as follows: SCD =1 2(3<cos2θCD >−1),(1) where θCD is the angle between the direction normal to the surface of the membrane and a CH-bond. It is remarkable that SCD can be obtained from 2H NMR experiments [46], making it a suitable property to verify the reliability of the simulations. The averaged results are shown in Fig. 3 for both tail chains of several classes of lipids considered in this work. The results indicate that all classes of lipids investigated show maxima of SCD between 0.37 and 0.43 units, namely with lower values found at 323 K and higher at 310 K. These values larger than 0.4 are in overall good agreement with the results of Ing´ olfsson et al. [31] for the full brain cell CG model reported by such authors. Since SCD is an indication of the degree of ordering (linearity) of the acyl chains, we confirm that the acyl chains of all types of lipids exhibit slightly greater ordering at 310 K, which becomes more disordered at 323 K, as expected. On the other hand, we can observe two additional features: (1) the ordering at Fig. 2. Backbone structures. Sketches of molecular structures of TMZ and the main classes of lipids considered in this work: cholesterol (CHOL), di-palmytoilphosphatidylcholine (DPPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE), 1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine (POPS), N-(hexadecanoyl)-hexadecasphing-4-enine-1-phosphocholine (DPSM) and 1,2-dipalmitoyl-sn-glycero-3-succinate (DPGS). The latter are cerebrosides and they are lipids exclusively located at the brain cells. Part of hydrogen‑carbon bonds are not shown, for the sake of clarity. The highlighted sites will be referred to in the text by the same labels. Table 1 Composition of the two leaflets (outer, inner) of the brain membrane, with percentages and number of lipids N. We should note that DPGS lipids are only located at the outer layer, whereas POPS lipids are only located in the inner layer of the BPM. Class Species Outer layer (%) Nouter Inner layer (%) Ninner PC DPPC 24.2 48 13.6 27 PE POPE 11.0 22 21.3 43 SM DPSM 11.0 22 6.0 12 PS POPS 0.0 0 14.5 29 cerebrosides DPGS 9.4 19 0.0 0 CHOL CHOL 44.4 89 44.6 89 Total –100.0 200 100.0 200 Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 3 the inner layer (triangles) is slightly greater than at the outer layer (circles) which indicates a more rigid phase for the inner leaflet and (2) those lipids having unsaturated bonds in one of their acyl tails, such as POPE, DPSM, POPS and DPSG show highly disordered profiles of SCD at the double bond positions. Other relevant parameters for the control of the simulations are the area per lipid A and the thickness Δz of the membrane. A is usually obtained considering the surface of the membrane along the XY plane divided by the number of lipids and cholesterol [47]. In the present case, the two leaflets contain the same number of lipids and their XY plane surfaces are the same at each instant of the simulation, but since the classes of lipids in each membrane are different (see Table 1), small differences could arise. As a rough approach, the averaged areas per lipid of two membranes (at lowest and highest temperatures) were computed geometrically as described above and are reported in Table 2. We obtained values of A around 41 Å2 at 310 K and slightly higher at 323 K, indicating a larger thermal disorder in the latter case. The full profile of A as a function of the simulation time is reported in Fig. S1 of SD. Our results are in good overall agreement with other computational works, such as the reference work of Ing´ olfsson et al. [31] or the study of Yee et al. [48] on the BPM of a healthy brain compared to another having Alzheimer’s disease. The values of A by these authors were 46 and 42.3Å2 respectively, i.e. with discrepancies between 3 and 11 %, probably due to the different models considered. Nevertheless, in order to obtain significantly more accurate values, areas per lipid can be obtained using the method of Voronoi tessellation, which is able to compute the contribution of each class of lipid to the averaged value of A. For full details of the calculations, see Refs. [49–50]. In the present work, we have employed the script “scipy. spatial.Voronoi” [51] combined with the tool “Convex Hull” [52]. The computed values for each leaflet and at the two extreme temperatures are reported in Table2. There we can notice that the overall, weighted average values are in excellent agreement with the previous geometrical Fig. 3. Order parameter. ∣SCD∣ for two acyl tails (sn1 and sn2 shown in Fig. 2) of selected lipid species considered in the membrane model of brain cells: DPPC (A-B); POPE (C – D); DPSM (E-F); POPS (G) and DPGS (H). Lipids at the inner layer are indicated with triangles and those at the outer layer are represented with circles. POPS in panel G is present only in the inner layer, while DPGS in panel H is only in the outer layer, and the circle represents sn1 while the triangle represents sn2. Systems are also represented at temperatures of 310 K (red) and 323 K (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 4 calculation, but using this methodology, we can separately extract the contribution of each class of lipid. In this way, we should notice that: (1) the overall effect of rising temperature is a moderate growth of the values of A for all lipid types; (2) lipids of the same class can contribute with slightly different areas per lipid depending of their surroundings, which are different in the two leaflets; (3) the effect of cholesterol (44–45 % of each leaflet) is totally dominant over the averaged values. In passing, we can observe that the numbers of the areas per lipid corresponding to specific species reported by Greenfield et al. [50] for phospholipids of the phosphatidylglycerol and lysylphosphatidylglycerol classes are quite similar to those reported here, i. e. around 40Å2. The thickness of the membrane may provide additional clues about the influence of cholesterol on the mechanical properties of plasma membranes, such as rigidity and capability of allowing the movement of species in and out of the cell. We have obtained the thickness of the membrane Δz (defined and reported as a function of time in Fig. S2 of SD). The results reported here (around 42 Åat 310 K) are again qualitatively close to those reported in Refs. [31] (41 Å) and Ref. [48] (47.3 Å) where more detailed models were considered, but only reported at 310 K. The full profile of Δz as a function of the simulation time is reported in Fig. S2 of SD, showing fluctuations of the order of 5 % around the averaged values, reported in Table2, where we can observe that as the temperature increases from 310 K to 323 K, thickness slightly decreases, while area per lipid increases. These results indicate a slightly higher degree of disorder as temperature increases, given the higher values of A and smaller values of Δz at 323 K as expected. 3.2. Radial distribution functions of TMZ bound to cholesterol and selected lipids A direct route to the characterisation of the local structure of each atomic species of the system is usually obtained by means of normalised radial distribution functions (RDF) gAB(r)for two different species A and B (see Eq.(2) of Ref. [45], for instance), indicating the probability of finding a species A at a given distance r of a species B. Among the wide variety of possible RDF that could be computed, we have considered only the most relevant RDF based on the first coordination shells of TMZ with cholesterol and selected lipids. The remaining RDF indicate low maxima at distances significantly longer than the typical hydrogen bond (HB) values or show rather noisy profiles which indicate that the corresponding local structures are not stable enough. The selected RDF are reported in Fig. 4. As a general observation, we note a well-defined first coordination shell with a sharp peak located around 1.8–2.0 Å, which is indicative of HB between TMZ and the surrounding species. This range of distances is typical for HB interactions between small molecules and organic species (see Ref. [53] and references therein). Among the various lipids, the first peak in the radial distribution function (RDF) for the interaction Table 2 Area per lipid (in Å2) and membrane thickness Δz (in Å) at 310 K and 323 K. Estimated errors, computed from standard deviations, are in parentheses. Each lipid class has been computed using Voronoi tessellation. Property 310 K 323 K Thickness (Å) 42.6 (0.2) 42.1 (0.2) Area per lipid (Å2)41.2 (0.2) 42.4 (0.3) Area per lipid (in Å2, outer leaflet) DPGS 39.1 (0.7) 41.7 (1.2) DPPC 47.5 (0.2) 47.8 (0.1) POPE 45.8 (0.6) 46.8 (0.3) DPSM 40.5 (0.4) 42.0 (0.5) CHOL 37.4 (0.2) 38.7 (0.1) Weighted average 41.3 (0.1) 42.4 (0.1) Area per lipid (in Å2, inner leaflet) POPS 43.5 (1.3) 45.0 (1.2) DPPC 44.9 (0.1) 46.7 (0.2) POPE 44.6 (0.2) 46.1 (0.3) DPSM 42.7 (0.7) 42.2 (0.4) CHOL 37.6 (0.4) 38.7 (0.3) Weigthed average 41.2 (0.1) 42.5 (0.1) Fig. 4. Radial distribution functions. RDF for different atoms in TMZ with several lipids at 310 K (black) and 323 K (red) systems. In panel A, gN5tmz−Hchol(r)refers to the RDF of the atom named N5 (see Fig. 2) in TMZ and the hydrogen of the hydroxyl group of cholesterol. In panels B – F the selected atoms are also indicated with the same labels as in Fig. 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 5 between N5 of TMZ and the POPE/POPS lipids is notably higher than for the others, especially for the hydrogens at the headgroups of POPE/ POPS (see Fig. 4, panels C and D). This suggests stronger binding of TMZ to phosphatidylethanolamine (POPE) and phosphatidylserine (POPS) lipids. Conversely, panels B and E in Fig. 4 show that DPGS, a specific cerebroside characteristic of the outer leaflet of the brain cell membrane, also contributes to TMZ binding through hydrogens H3 and H4. Notably, a relatively stronger interaction between TMZ and DPGS is observed at the higher temperature (323K). The noisy profiles in panels B and E indicate that the TMZ-DPGS interaction has a shorter lifetime compared to the TMZ-POPE/POPS binding, as long-lived interactions (such as HB) result in smoother RDF profiles, whereas short-lived interactions undergo frequent breaking and reforming events. In the particular case of TMZ-DPGS hydrogen-bonding we can observe mainly HB formed by nitrogen N5 of TMZ with H3 and H4 from DPGS and an oxygen‑hydrogen interaction of H1 of TMZ with O2 in POPS, but only well defined at 310 K. Comparing systems at 310 K and 323 K, it can be noticed that temperature plays a minor role in most cases. In a few cases, such as TMZ-DGPS bindings (see panels B and E of Fig. 4) we can notice an enhancement of the first maxima at 323 K, which can be attributed to more frequent and durable events of TMZDPGS pairings than at 310 K, when TMZ shows a bigger affinity to bind cholesterol (panel A in Fig. 4) and POPS lipids (panel F of Fig. 4). 3.3. Estimation of free energy barriers for TMZ crossing From a general perspective, the calculation of the Helmholtz or Gibbs free energy differences for binding processes or configurational changes is a difficult task and it requires a considerable amount of computer time and a precise knowledge of the hypersurface of the potential energy of the system [54]. This can be explored by means of methods such as metadynamics [55–57], hybrid quantum mechanics/molecular mechanics methods [58] or transition path sampling [59–62]. However, a usual way to obtain free energy estimations is through the calculation of potentials of mean force (PMF) through a variety of methods such as umbrella sampling [63], constrained molecular dynamics [64] or adaptive biasing force methods, among many others. In this Section we computed free energy profiles by means of the adaptive biasing force (ABF) method [65,66], where the PMF is derived from integrating the mean force Fz over the reaction coordinates ξ0 [67–69]. The total work through the membrane is: ΔW= − [∫ξfinal ξoriginal 〈Fz(ξ0) 〉dξ0],(2) where ξoriginal is the original reaction coordinate of the TMZ and ξfinal is the final reaction coordinate of the TMZ. 〈Fz(ξ0) 〉 is the average force in the z-direction exerted on the TMZ at ξ0. The ABF method offers the advantage of not requiring the estimation of free energy barriers and providing more diverse sampling. In typical simulations, the substrate is driven in a specific direction by the total force acting on it. With ABF, a biasing force in the opposite direction is introduced, creating a balance of the substrate and leading to random moving. Nevertheless, ABF methods may exhibit inaccuracies arising from insufficient sampling and sensitivity to the choice of reaction coordinate, which can hinder convergence in complex systems [68,70,71]. The reaction coordinates in the present case represent the distance along the z-axis from the center of mass of TMZ to a “dummy atom” located at the point (0.0, 0.0, 0.0). The reaction coordinates are determined by introducing a dummy atom instead of real molecules or atoms to prevent deviation caused by their movement. These collective variables in small size improve the performance of ABF calculation. The reaction coordinates are divided into several windows. TMZ will be confined within specific windows by harmonic forces in the Z-axis direction at the boundaries, which enhances the sampling efficiency. Three TMZ membrane-penetrating ABF systems were set up at 310 K, 323 K and an additional intermediate temperature at 316.5 K, to gain a comprehensive understanding of the impact of system temperature. Each window has a width of 5Å, and the PMF is calculated every 0.1Å. The force constant at the window boundaries is 20 kcal/mol/Å2. A sample count of the times that TMZ has been placed at different positions in the system in one ABF calculation for each temperature is reported in Fig. S3 of SD, in order to show the reliability of the ABF calculations and that sufficent sampling was achieved in all cases. The averaged results for the two free energy barriers are shown in Fig. 5 and the four component individual ABF trajectories used to obtain the averages and error bars are reported in Fig.S4 of SD. As general features, we can observe that: 1. All PMF profiles show a free energy maximum when TMZ reaches the center of the membrane, regardless of the direction taken either starting from the external aqueous solution or from the internal cytoplasm; 2. Fluctuations of all free energy profiles are quite large at all temperatures, especially when TMZ is close to the center of the membrane, indicating important configurational changes, mainly related to the different interactions of TMZ with surrounding lipids; 3. At the lowest temperature (310 K, i.e. normal human body temperature) the free energy barrier required by TMZ to reach the center of the BPM is the largest, of the order of 50 ±5 kcal/mol; 4. The two sets of free energy profiles (i.e. when TMZ is crossing outer and inner layers) are very similar each other with an almost perfect match between the two directions, but showing a small bias to larger values (of around 1–3 kcal/mol) when TMZ is moved from the internal solution and forced to cross the inner layer of the BPM; 5. We observe a clear monotonic behaviour, with larger free energy barriers at the lowest temperature and a tendency to decrease when temperature is raised. In addition, we have found a neat reduction of a factor ∼1.8–1.4 (from 49 to 28 kcal/mol to cross the outer layer and from 50 to 35 to cross the inner layer) when the temperature of the system is raised from 310 to 323 K. This would indicate a significantly better ability of TMZ to cross the BPM when temperature is higher than the normal one; 6. As the temperature rises, the difference in the free energy barrier that TMZ needs to surpass to enter/leave the BPM increases: for instance, at 323 K, the barrier grows from 28 kcal/mol to cross the outer layer to 35 kcal/mol to cross the inner layer, making it progressively more difficult for TMZ to be expelled from the BPM once it has entered. We can rationalize these findings as follows. First, regarding the large error bars observed in the ABF simulations, we explored the preferential conformations of TMZ in the case of the largest error bars (310K), where TMZ moves from the external solution to the center of the BPM through the outer layer, as seen in the blue line of Fig. 5A. The main interactions of TMZ for the two extreme ABF profiles are reported in Fig. S5 of the Supporting Information (SD). In the highest free energy profile (blue line), TMZ forms three strong hydrogen bonds (HB) with two cholesterol molecules, whereas in the lowest free energy profile (red line), TMZ forms two HBs with two DPPC molecules. This suggests that the strongest bindings of TMZ can vary substantially across different ABF trajectories, resulting in quantitatively different free energy barriers. On the other hand, there is also the possibility of TMZ undergoing recrossings, i.e., moving backward rather than crossing the full bilayer, returning to the solution where the trajectory started. The small differences in the free energy barriers reported in Fig. 5 suggest that both forward and backward motions starting near the center of the BPM are equally probable. To investigate this, we ran four independent unbiased MD simulations where TMZ was initially placed near the center of the BPM (z∼0 Å), using the coordinates from the ABF trajectories. The results for TMZ’s position along the z-axis are shown in Fig. S6 of the SD. Although the statistics for these auxiliary MD trajectories are limited, the goal was not to compute precise crossing rates (which would require Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 6 running hundreds of independent simulations) but to assess the possibility of recrossings. In two of the four trajectories, TMZ crossed the outer layer towards the “external” solution within 20 ns. In the remaining two trajectories, TMZ crossed the inner layer, taking between 40 and 50 ns. This suggests a stronger interaction between TMZ and the headgroups of lipids in the inner leaflet than with those in the outer leaflet. These findings are consistent with the results shown in Fig. S4 of the SD, where it is slightly more energetically costly for TMZ to cross the inner leaflet interface than the outer leaflet. 4. Summary and conclusions In the present study, we conducted a comprehensive computational analysis of TMZ permeation across a brain-like plasma membrane bilayer. We introduced a TMZ molecule at the interface of a model BPM and systematically observed its interactions with solvent water, ions, and the distinct regions of the membrane, which include both the outer and inner layers. Our results confirm the proper equilibration of the system, and structural data indicate that the characteristic order parameter used to verify the system’s liquid-ordered state, as well as other properties such as the area per lipid and membrane thickness, show good qualitative agreement with results reported by Ing´ olfsson et al. [31] and Yee et al. [48], where slightly different BPM models were employed. In addition, we have observed that the model BPM membrane designed in the present work has a nearly planar geometry for both leaflets, without any significantly large curvatures or stress on the membrane (see Fig. S7 of SD) observed during the whole time range (0,5 μ s) of the simulations. Our free energy profiles for TMZ transmembrane crossing obtained from adaptive biasing force methods indicate that TMZ has the highest affinity for staying at the interface between the BPM and the aqueous solution, rather than spontaneously crossing the membrane via simple diffusion. This suggests that in the real system, some agents such as membrane ABC transporters [72] or SLC transporters [73,74] may play a key role in the crossing of TMZ through the external layer of the brain cell membrane. Interestingly, when TMZ has reached the internal hydrophobic regions of the BPM, small fluctuations in its position can lead to: (1) a successful crossing to the internal regions and eventually reach the first cytoplasmic region of the brain or (2) to recrossings towards the external interface of the BPM and eventual return to the external aqueous solution. Our ABF calculations indicate that TMZ needs about 50 kcal/mol to cross the BPM at the normal human body temperature of 310 K, whereas this barrier decreases to about 28 kcal/mol at 323 K. Data obtained by RDF indicate that the formation of strong TMZcholesterol, TMZ-POPS, TMZ-POPE and TMZ-DPGS pairings are the main responsible for TMZ hydrogen bonding interactions that drive the permeation of TMZ through the BPM barrier. Finally, we would like to remark that optimizing the molecular design of TMZ through advanced drug delivery systems can enhance its pharmacokinetic properties. This is normally done by improving stability, prolonging circulation time, and increasing brain bioavailability. For instance, TMZ-loaded PEGylated liposomes have demonstrated a 4.2-fold increase in brain drug concentration and a 1.6-fold increase in blood AUC compared to the free drug [75]. Similarly, TMZ encapsulated in niosomal formulations has shown enhanced brain targeting and reduced systemic exposure, leading to improved therapeutic efficacy [76]. As an alternative, from our results we observe that a reduction of the size of the free energy barrier that TMZ needs to cross to access the brain could be carried out by reducing its ability to form HB with lipids. This effect could be achieved by replacing or masking HB donor or acceptor groups (for instance by converting hydroxyl or amide groups into less polar functionalities or by introducing steric hindrance around HB-prone sites. In summary, these changes can decrease nonspecific interactions with membrane lipids, potentially improving membrane permeability and pharmacokinetics. CRediT authorship contribution statement Yanhong Ge: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. Huixia Lu: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Methodology, Investigation, Formal analysis, Conceptualization. Jordi Martí: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Methodology, Investigation, Formal analysis, Conceptualization. Fig. 5. Free energy profiles of TMZ crossing the brain cell membrane bilayer. (A): TMZ from solution towards the center of the BPM, crossing outer layer, (B): TMZ from the solution towards the center of the BPM, crossing the inner layer. The origin of the free energy is taken for TMZ being solvated in the aqueous solution, for each of the three temperatures. Each PMF curve is the average result obtained from four independent ABF simulations, so that the shadowed areas indicate error bars and they are obtained as the largest deviations between the average and the component profiles. Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 7 Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors acknowledge financial support provided by the Spanish Ministry of Science, Innovation and Universities. Huixia Lu thanks the financial support by the “Margarita Salas” grant which is funded by the European Union–NextGenerationEU. This publication is a part of the I + D +i project with reference PID2021-124297NB-C32, founded by MCIN/AEI/10.13-039/501100011033 and “FEDER Una manera de hacer Europa”. Yanhong Ge is a Ph.D. fellow from the China Scholarship Council (grant 202306230043). J.M. thanks the Generalitat de Catalunya for the support through the grant Grup de Recerca SGR-Cat2021 Condensed, Complex and Quantum Matter Group reference 2021SGR01411 and to the Polytechnic University of Catalonia-Barcelona Tech through the funding AGRUPS. The authors thankfully acknowledge computer resources from MareNostrum5 supercomputer as well as technical support provided by BSC (RES-BCV-2024-2-0006). Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.bpc.2025.107457. Data availability All related files including an example used to calculate the binding free energy of this work can be found in our repository of Zenodo: https://doi.org/10.5281/zenodo.12745251. References [1] H.S. Friedman, T. Kerby, H. Calvert, Temozolomide and treatment of malignant glioma, Clin. Cancer Res. 6 (7) (2000) 2585–2597. [2] A.C. Tan, D.M. Ashley, G.Y. L´ opez, M. Malinzak, H.S. Friedman, M. Khasraw, Management of glioblastoma: state of the art and future directions, CA Cancer J. Clin. 70 (4) (2020) 299–312. [3] W.A. Yung, M.D. Prados, R. Yaya-Tur, S.S. Rosenfeld, M. Brada, H.S. Friedman, R. Albright, J. Olson, S.M. Chang, A.M. O’Neill, et al., Multicenter phase ii trial of temozolomide in patients with anaplastic astrocytoma or anaplastic oligoastrocytoma at first relapse, J. Clin. Oncol. 17 (9) (1999) 2762–2771. [4] A. Thomas, M. Tanaka, J. Trepel, W.C. Reinhold, V.N. Rajapakse, Y. Pommier, Temozolomide in the era of precision medicine, Cancer Res. 77 (4) (2017) 823–826. [5] D. Wu, Q. Chen, X. Chen, F. Han, Z. Chen, Y. Wang, The blood–brain barrier: structure, regulation, and drug delivery, Signal Transduct. Target. Ther. 8 (1) (2023) 217. [6] W. Risau, H. Wolburg, Development of the blood-brain barrier, Trends Neurosci. 13 (5) (1990) 174–178. [7] P. Ballabh, A. Braun, M. Nedergaard, The blood–brain barrier: an overview: structure, regulation, and clinical implications, Neurobiol. Dis. 16 (1) (2004) 1–13. [8] O. Van Tellingen, B. Yetkin-Arik, M. De Gooijer, P. Wesseling, T. Wurdinger, H. De Vries, Overcoming the blood–brain tumor barrier for effective glioblastoma treatment, Drug Resist. Updat. 19 (2015) 1–12. [9] M. Rubio-Camacho, J.A. Encinar, M.J. Martínez-Tom´ e, R. Esquembre, C.R. Mateo, The interaction of temozolomide with blood components suggests the potential use of human serum albumin as a biomimetic carrier for the drug, Biomolecules 10 (7) (2020) 1015. [10] H. Strobel, T. Baisch, R. Fitzel, K. Schilberg, M.D. Siegelin, G. Karpel-Massler, K.- M. Debatin, M.-A. Westhoff, Temozolomide and other alkylating agents in glioblastoma therapy, Biomedicines 7 (3) (2019) 69. [11] D. Petrenko, V. Chubarev, N. Syzrantsev, N. Ismail, V. Merkulov, S. Sologova, E. Grigorevskikh, E. Smolyarchuk, R. Alyautdin, Temozolomide efficacy and metabolism: the implicit relevance of nanoscale delivery systems, Molecules 27 (11) (2022) 3507. [12] I.C. Lopes, S.C.B. de Oliveira, A.M. Oliveira-Brett, Temozolomide chemical degradation to 5-aminoimidazole-4-carboxamide–electrochemical study, J. Electroanal. Chem. 704 (2013) 183–189. [13] M.J. Tisdale, Antitumour imidazotetrazines—xv: role of guanine o6 alkylation in the mechanism of cytotoxicity of imidazotetrazinones, Biochem. Pharmacol. 36 (4) (1987) 457–462. [14] B.J. Denny, R.T. Wheelhouse, M.F. Stevens, L.L. Tsang, J.A. Slack, Nmr and molecular modeling investigation of the mechanism of activation of the antitumor drug temozolomide and its interaction with dna, Biochemistry 33 (31) (1994) 9045–9051. [15] J. Zhang, M.F.G. Stevens, T.D. Bradshaw, Temozolomide: mechanisms of action, repair and resistance, Curr. Mol. Pharmacol. 5 (1) (2012) 102–114. [16] K. Kelly, M. Sacapano, K. Prosolovich, J. Ong, K. Black, Blood brain barrier permeability to temozolomide, Cancer Res. 65 (9 Supplement) (2005) 330. [17] M.C. de Gooijer, N.A. de Vries, T. Buckle, L.C. Buil, J.H. Beijnen, W. Boogerd, O. van Tellingen, Improved brain penetration and antitumor efficacy of temozolomide by inhibition of abcb1 and abcg2, Neoplasia 20 (7) (2018) 710–720. [18] G. Blackledge, J. Roberts, S. Kaye, R. Taylor, J. Williams, B. De Stavola, B. Uscinska, A phase ii study of mitozolomide in metastatic transitional cell carcinoma of the bladder, Eur. J. Cancer Clin. Oncol. 25 (2) (1989) 391–392. [19] K. Bouzinab, H.S. Summers, M.F. Stevens, C.J. Moody, N.R. Thomas, P. Gershkovich, N. Weston, M.B. Ashford, T.D. Bradshaw, L. Turyanska, Delivery of temozolomide and n3-propargyl analog to brain tumors using an apoferritin nanocage, ACS Appl. Mater. Interfaces 12 (11) (2020) 12609–12617. [20] M. Chakravarty, P. Ganguli, M. Murahari, R.R. Sarkar, G.J. Peters, Y. Mayur, Study of combinatorial drug synergy of novel acridone derivatives with temozolomide using in-silico and in-vitro methods in the treatment of drug-resistant glioma, Front. Oncol. 11 (2021) 625899. [21] S.D. Precilla, S.S. Kuduvalli, E.A. Praveena, S. Thangavel, T. Anitha, Integration of synthetic and natural derivatives revives the therapeutic potential of temozolomide against glioma-an in vitro and in vivo perspective, Life Sci. 301 (2022) 120609. [22] R. Li, D. Tang, J. Zhang, J. Wu, L. Wang, J. Dong, The temozolomide derivative 2tp400 inhibits glioma growth via administration route of intravenous injection, J. Neuro-Oncol. 116 (2014) 25–30. [23] J. Yin, X. Wang, X. Ge, F. Ding, Z. Shi, Z. Ge, G. Huang, N. Zhao, D. Chen, J. Zhang, et al., Hypoxanthine phosphoribosyl transferase 1 metabolizes temozolomide to activate ampk for driving chemoresistance of glioblastomas, Nat. Commun. 14 (1) (2023) 5913. [24] M.J. Ramalho, S. Andrade, M.´ A.N. Coelho, J.A. Loureiro, M.C. Pereira, Biophysical interaction of temozolomide and its active metabolite with biomembrane models: the relevance of drug-membrane interaction for glioblastoma multiforme therapy, Eur. J. Pharm. Biopharm. 136 (2019) 156–163. [25] I. Kracun, H. Rosner, V. Drnovsek, M. Heffer-Lauc, C. Cosovi´ c, G. Lauc, Human brain gangliosides in development, aging and disease, Int. J. Dev. Biol. 35 (3) (1991) 289–295. [26] L. Svennerholm, K. Bostr¨ om, B. Jungbjer, L. Olsson, Membrane lipids of adult human brain: lipid composition of frontal and temporal lobe in subjects of age 20 to 100 years, J. Neurochem. 63 (5) (1994) 1802–1811. [27] M. Aureli, S. Grassi, S. Prioni, S. Sonnino, A. Prinetti, Lipid membrane domains in the brain, Biochim. et Biophys. Acta (BBA)-Mol. Cell Biolo. Lipids 1851 (8) (2015) 1006–1016. [28] Y. Arai, J.L. Sampaio, M. Wilsch-Br¨ auninger, A.W. Ettinger, C. Haffner, W. B. Huttner, Lipidome of midbody released from neural stem and progenitor cells during mammalian cortical neurogenesis, Front. Cell. Neurosci. 9 (2015) 325. [29] P. Sj¨ ovall, J. Lausmaa, B. Johansson, Mass spectrometric imaging of lipids in brain tissue, Anal. Chem. 76 (15) (2004) 4271–4278. [30] I. Chakrabarti, M. Cockburn, W. Cozen, Y.-P. Wang, S. Preston-Martin, A population-based description of glioblastoma multiforme in los Angeles county, 1974–1999, Cancer, Interdisciplinary International Journal of the American Cancer Society 104 (12) (2005) 2798–2806. [31] H.I. Ing´ olfsson, T.S. Carpenter, H. Bhatia, P.-T. Bremer, S.J. Marrink, F. C. Lightstone, Computational lipidomics of the neuronal plasma membrane, Biophys. J. 113 (10) (2017) 2271–2280. [32] S. Jo, T. Kim, V.G. Iyer, W. Im, Charmm-gui: a web-based graphical user interface for charmm, J. Comput. Chem. 29 (11) (2008) 1859–1865. [33] S. Jo, J.B. Lim, J.B. Klauda, W. Im, Charmm-gui membrane builder for mixed bilayers and its application to yeast membranes, Biophys. J. 97 (1) (2009) 50–58. [34] J. Yang, J. Martí, C. Calero, Pair interactions among ternary dppc/popc/cholesterol mixtures in liquid-ordered and liquid-disordered phases, Soft Matter 12 (20) (2016) 4557–4561. [35] J.C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with namd, J. Comput. Chem. 26 (16) (2005) 1781–1802. [36] J.B. Klauda, R.M. Venable, J.A. Freites, J.W. O’Connor, D.J. Tobias, C. MondragonRamirez, I. Vorobyov, A.D. MacKerell Jr., R.W. Pastor, Update of the charmm allatom additive force field for lipids: validation on six lipid types, J. Phys. Chem. B 114 (23) (2010) 7830–7843. [37] J.B. Lim, B. Rogaski, J.B. Klauda, Update of the cholesterol force field parameters in charmm, J. Phys. Chem. B 116 (1) (2012) 203–210. [38] J.A. Lemkul, J. Huang, B. Roux, A.D. MacKerell Jr., An empirical polarizable force field based on the classical drude oscillator model: development history and recent applications, Chem. Rev. 116 (9) (2016) 4983–5013. [39] J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B.L. De Groot, H. Grubmüller, A.D. MacKerell Jr., Charmm36m: an improved force field for folded and intrinsically disordered proteins, Nat. Methods 14 (1) (2017) 71–73. [40] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, A smooth particle mesh ewald method, J. Chem. Phys. 103 (19) (1995) 8577–8593. [41] H.J. Berendsen, J.V. Postma, W.F. Van Gunsteren, A. Di Nola, J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (8) (1984) 3684–3690. [42] G.J. Martyna, D.J. Tobias, M.L. Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys. 101 (5) (1994) 4177–4189. Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 8 [43] S.E. Feller, Y. Zhang, R.W. Pastor, B.R. Brooks, Constant pressure molecular dynamics simulation: the langevin piston method, J. Chem. Phys. 103 (11) (1995) 4613–4621. [44] G.W. Stockton, I.C. Smith, A deuterium nuclear magnetic resonance study of the condensing effect of cholesterol on egg phosphatidylcholine bilayer membranes. i. Perdeuterated fatty acid probes, Chem. Phys. Lipids 17 (2–3) (1976) 251–263. [45] H. Lu, J. Martí, Effects of cholesterol on the binding of the precursor neurotransmitter tryptophan to zwitterionic membranes, J. Chem. Phys. 149 (16) (2025). [46] P.L. Yeagle, The Membranes of Cells, Academic Press, 2016. [47] P.R. Pandey, S. Roy, Headgroup mediated water insertion into the dppc bilayer: a molecular dynamics study, J. Phys. Chem. B 115 (12) (2011) 3155–3163. [48] S.M. Yee, R.J. Gillams, S.E. McLain, C.D. Lorenz, Effects of lipid heterogeneity on model human brain lipid membranes, Soft Matter 17 (1) (2021) 126–135. [49] W. Shinoda, S. Okazaki, A voronoi analysis of lipid area fluctuation in a bilayer, J. Chem. Phys. 109 (4) (1998) 1517–1521. [50] M.L. Greenfield, L.M. Martin, F. Joodaki, Computing individual area per head group reveals lipid bilayer dynamics, J. Phys. Chem. B 126 (50) (2022) 10697–10711. [51] H. Yin, G. Cusatis, Ringspy: a python package for voronoi mesh generation of cellular solids with radial growth pattern, J. Open Source Softw. 8 (83) (2023) 4945. [52] E.W. Weisstein, Convex hull. https://mathworld.wolfram.com/, 2025. [53] J. Martí, H. Lu, Microscopic interactions of melatonin, serotonin and tryptophan with zwitterionic phospholipid membranes, Int. J. Mol. Sci. 22 (6) (2021) 2842. [54] F. M. Ytreberg, R. H. Swendsen, D. M. Zuckerman, Comparison of free energy methods for molecular systems, J. Chem. Phys. 125 (18). [55] A. Laio, M. Parrinello, Escaping free-energy minima, Proc. Natl. Acad. Sci. 99 (20) (2002) 12562–12566. [56] J. Yang, C. Calero, M. Bonomi, J. Martí, Specific ion binding at phospholipid membrane surfaces, J. Chem. Theory Comput. 11 (9) (2015) 4495–4499. [57] H. Lu, J. Marti, Cellular absorption of small molecules: free energy landscapes of melatonin binding at phospholipid membranes, Sci. Rep. 10 (1) (2020) 9235. [58] H.M. Senn, W. Thiel, Qm/mm methods for biological systems, Atomist. Approach. Mode. Biol.: Quant. Chem. to Mol. Simulat. (2007) 173–290. [59] J. Martí, F.S. Csajka, D. Chandler, Stochastic transition pathways in the aqueous sodium chloride dissociation process, Chem. Phys. Lett. 328 (1–2) (2000) 169–176. [60] P.G. Bolhuis, D. Chandler, C. Dellago, P.L. Geissler, Transition path sampling: throwing ropes over rough mountain passes, in the dark, Annu. Rev. Phys. Chem. 53 (1) (2002) 291–318. [61] J. Martí, F.S. Csajka, Transition path sampling study of flip-flop transitions in model lipid bilayer membranes, Phys. Rev. E 69 (6) (2004) 061918. [62] C. Dellago, P.G. Bolhuis, Transition path sampling and other advanced simulation techniques for rare events, in: Advanced computer simulation approaches for soft matter sciences III, 2009, pp. 167–233. [63] J. K¨ astner, Umbrella sampling, Wiley interdisciplinary reviews: computational molecular, Science 1 (6) (2011) 932–942. [64] M. Sprik, G. Ciccotti, Free energy from constrained molecular dynamics, J. Chem. Phys. 109 (18) (1998) 7737–7744. [65] J. H´ enin, J. Gumbart, C. Chipot, Free Energy Calculations Along a Reaction Coordinate: A Tutorial for Adaptive Biasing Force Simulations, 2025. [66] J. H´ enin, C. Chipot, Overcoming free energy barriers using unconstrained molecular dynamics simulations, J. Chem. Phys. 121 (7) (2004) 2904–2914. [67] E. Darve, D. Rodríguez-G´ omez, A. Pohorille, Adaptive biasing force method for scalar and vector free energy calculations, J. Chem. Phys. 128 (14) (2025). [68] E. Darve, A. Pohorille, Calculating free energies using average force, J. Chem. Phys. 115 (20) (2001) 9169–9183. [69] C. Chipot, J. H´ enin, Exploring the free-energy landscape of a short peptide using an average force, J. Chem. Phys. 123 (24) (2025). [70] C. Chipot, A. Pohorille, Free Energy Calculations: Theory and Applications in Chemistry and Biology, 86 of Springer Series in Chemical Physics, Springer, 2007. [71] J. Comer, J.C. Gumbart, J. H´ enin, T. Leli` evre, A. Pohorille, C. Chipot, The adaptive biasing force method: everything you always wanted to know but were afraid to ask, J. Phys. Chem. B 119 (3) (2015) 1129–1151. [72] D.C. Rees, E. Johnson, O. Lewinson, Abc transporters: the power to change, Nat. Rev. Mol. Cell Biol. 10 (3) (2009) 218–227. [73] L. Lin, S.W. Yee, R.B. Kim, K.M. Giacomini, Slc transporters as therapeutic targets: emerging opportunities, Nat. Rev. Drug Discov. 14 (8) (2015) 543–560. [74] S.K. Nigam, What do drug transporters really do? Nat. Rev. Drug Discov. 14 (1) (2015) 29–44. [75] J. Vanza, P. Jani, N. Pandya, H. Tandel, Formulation and statistical optimization of intravenous temozolomide-loaded pegylated liposomes to treat glioblastoma multiforme by three-level factorial design, Drug Dev. Ind. Pharm. 44 (6) (2018) 923–933. [76] A. De, N. Venkatesh, M. Senthil, B.K.R. Sanapalli, R. Shanmugham, V.V.S.R. Karri, Smart niosomes of temozolomide for enhancement of brain targeting, Nanobiomedicine 5 (2018), 1849543518805355. Y. Ge et al. Biophysical Chemistry 324 (2025) 107457 9