Membrane Permeability in Biological Systems : A Systems Biology Perspective

Domains of modern computational science are converging on complex problems in the general field of systems biology. There is now a credible possibility of modeling drug delivery vesicles (liposomes) and their properties with qualitative and quantitative insight coming from atomistic calculations. Such work leads to extend accurate quantum mechanical (QM) methods to large-scale atomistic molecular mechanics (MM) modeling. Also, it allows coarse-grain parameterization and, further on, mathematical models developing in systems biology. Starting from Density Functional Theory (DFT) methods and recent theoretical studies of the properties of phospholipid molecules, consistent QM/MM and MM methodologies, as well as their counterparts as functions of time, BornOppenheimer dynamics (BOMD) and molecular dynamics (MD) can also be generated. These methods applied to understand and predict the properties of lipid layers (bilayers), in particular, fluidity, and their evolution in the presence of ions (physiological conditions) and of additional biological agents. The strategy is to build links between QM electronic-structure calculations at the finest molecular scale and equilibrium and nonequilibrium MM simulations to address mesoscopic system properties, using error control across the different scales. This strategy has not been used up to now in the domain of biochemical simulations, which explains the discrepancies found, for lipids, between atomistic QM and MM approaches. The reason for this situation comes from the fact that most of the MM parameterizations have been established on the basis of macroscopic experimental properties, most generally giving information at a time scale much larger than the simulation time of less than 50ns (for example, NMR time scale is between 10-3 to 10-4s according to the field frequency).The methodology incorporates in the MM approaches, as much as possible, the information provided by QM studies of the properties explored at the mesoscopic level. Thus, understanding the formation of vesicles from phospholipids bilayers and their fluidity and permeability properties is the basis of a large number of applications in the domain of drug delivery, in particular release of the active species according to the pH or ionic concentration changes. Prediction of structural changes (phase transition in particular) of membranes by modification of one or several constituents or addition of external molecular species may have potential therapeutic applications. Understanding the basic principles of biomembranes (lipid bilayers), which govern and mediate various biologically relevant processes on the microscopic level is one of the great challenges in biology. To investigate the characteristics of the membranes and to obtain the intriguing physicochemical aspects of membranes systems many experiments have been and are still being performed. Recent development of new algorithms and revolutionary advances in the computational power permits large-scale molecular dynamic (MD) simulations of interactions between biomembranes and small non-water polar molecules as well as simulations of two component mixtures of phospholipids membranes and other natural amphiphiles. Thus, allowing the development of a combination of computer simulations and experiments to analyze biomembrane properties with an even greater degree of detail encompasses our scientific endeavor.


Background
Domains of modern computational science are converging on complex problems in the general field of systems biology. There is now a credible possibility of modeling drug delivery vesicles (liposomes) and their properties with qualitative and quantitative insight coming from atomistic calculations. Such work leads to extend accurate quantum mechanical (QM) methods to large-scale atomistic molecular mechanics (MM) modeling. Also, it allows coarse-grain parameterization and, further on, mathematical models developing in systems biology. Starting from Density Functional Theory (DFT) methods and recent theoretical studies of the properties of phospholipid molecules, consistent QM/MM and MM methodologies, as well as their counterparts as functions of time, Born-Oppenheimer dynamics (BOMD) and molecular dynamics (MD) can also be generated. These methods applied to understand and predict the properties of lipid layers (bilayers), in particular, fluidity, and their evolution in the presence of ions (physiological conditions) and of additional biological agents. The strategy is to build links between QM electronic-structure calculations at the finest molecular scale and equilibrium and nonequilibrium MM simulations to address mesoscopic system properties, using error control across the different scales. This strategy has not been used up to now in the domain of biochemical simulations, which explains the discrepancies found, for lipids, between atomistic QM and MM approaches. The reason for this situation comes from the fact that most of the MM parameterizations have been established on the basis of macroscopic experimental properties, most generally giving information at a time scale much larger than the simulation time of less than 50ns (for example, NMR time scale is between 10 -3 to 10 -4 s according to the field frequency).The methodology incorporates in the MM approaches, as much as possible, the information provided by QM studies of the properties explored at the mesoscopic level. Thus, understanding the formation of vesicles from phospholipids bilayers and their fluidity and permeability properties is the basis of a large number of applications in the domain of drug delivery, in particular release of the active species according to the pH or ionic concentration changes. Prediction of structural changes (phase transition in particular) of membranes by modification of one or several constituents or addition of external molecular species may have potential therapeutic applications.
Understanding the basic principles of biomembranes (lipid bilayers), which govern and mediate various biologically relevant processes on the microscopic level is one of the great challenges in biology. To investigate the characteristics of the membranes and to obtain the intriguing physicochemical aspects of membranes systems many experiments have been and are still being performed. Recent development of new algorithms and revolutionary advances in the computational power permits large-scale molecular dynamic (MD) simulations of interactions between biomembranes and small non-water polar molecules as well as simulations of two component mixtures of phospholipids membranes and other natural amphiphiles. Thus, allowing the development of a combination of computer simulations and experiments to analyze biomembrane properties with an even greater degree of detail encompasses our scientific endeavor.
Molecular dynamic (MD) simulations are well suited for detailed analysis of the interactions between lipid bilayers and various small molecules, including water, chemicals, co-enzymes, peptides, oligonucleotides and proteins, as evidenced by the extensive body of published literature. Starting from Density Functional Theory (DFT) methods and recent theoretical studies of the properties of phospholipid molecules, consistent QM/MM and MM methodologies as well as their counterparts as function of time, Born-Oppenheimer dynamics (BOMD) and molecular dynamics (MD) could be explored. These methods will be applied to understand and predict the properties of lipid assemblies in water, in particular, fluidity, and their evolution in the presence of ions and biological agents. The main focus lies onto build links between QM electronic-structures and equilibrium and non-equilibrium MM simulations, using error control across the different scales. Moreover, there will be the integration of the methods and computer codes, particularly accounting for intermolecular polarization effects, which are of major importance in these systems.
Lipid bilayers have been studied extensively with classical MD since the early nineties. It is not possible to give an exhaustive review of these works, performed, most often, in conjunction with experimental results by electron microscopy, neutron diffraction, eventually single crystal X-Ray, NMR, IR spectroscopies, calorimetry, etc. Among the large number of MD studies devoted to phospholipid structural properties, dihedral angle values [1], head group flexibility [2], tail orientations [3], phase changes [4], hydration effects [5], interaction with ions [6] and molecules [7], evaluation of local order parameters [8,9] have been explored. Coarse grain methods have also been used to explore the formation and structure of small phospholipid vesicles [10].
Targeted delivery could bring a sharp improvement of existing approaches to diagnosis and treatment of various diseases. Therapeutic application of targeted drug delivery is limited though clinical trials of some immunotoxins and other antibody-drug conjugates are still being studied. Liposomal delivery systems offer major improvement in therapeutics through site specificity, their ability to escape from multi-drug resistance and the efficient delivery of targeted agent. Liposomes are potential drug carrier systems because of their ability to alter pharmacokinetics and reduced toxicity of drugs associated with them. In terms of targeting drug delivery by immunoliposomes, two anatomical targeting sites can be considered; divided into two phases: the transport phase in which the immunoliposomes travel from the site of administration to the target cells and the effector phase that includes the specific binding of immunoliposomes to the target cells with the subsequent delivery of entrapped drugs [11]. Henceforth, liposomes are correctly said as an "iceberg" to create novel individualized therapies. More recently, the introduction of liposomes serve as a viable delivery vehicle for chemotherapeutics, which include nanovectors, nanomaterials and microbubbles. One way of exploiting the diagnostic and therapeutic applications of microbubbles is with ultrasound. Ultrasound is one of the most common medical imaging methods used therapeutically [12]. Microbubbles comprising of spherical voids or cavities filled by a gas are stabilized by a coating material such as phospholipid, surfactant, denatured human serum albumin or synthetic polymer. Gas being less dense than liquids or solids, microbubbles comprises a pocket, region or structure of low density. This property of low density of microbubbles has a number of potentially important medical applications, including site-specific delivery, treatment of thrombosis and pulmonary delivery. Bioactive compounds can be incorporated into these carriers for site-specific delivery. Owing to the computational cost of quantum methods, energy calculations on large biological systems are mostly limited to the application of classical molecular mechanics. Some interactions however, such as those involving metal cofactors or the formation of covalent bonds, require the application of the laws of quantum mechanics.
Complexity of biology and the huge biological diversity among apparently similar individuals, mathematical models are clearly of fundamental importance for the effective development of tools to be used in personalized medicine. A number of researchers throughout the world have shown that ultrasound and microbubbles may have important applications in gene therapy [12] as there are different ways a microbubble can transport drugs. The mathematical modeling of biophysical phenomena is crucial for identifying the main parameters governing the spatio-temporal evolution of the system for elucidating their roles and quantifying their effects. Drugs may be attached to the membrane surrounding the microbubble or may be imbedded within the membrane itself. Materials like DNA may be bound noncovalently to the surface of the microbubbles; microbubbles might also be formulated to load the interior with drugs and gas or hydrophobic drugs can be incorporated into a layer of oily material that forms a film around the microbubble, which is then surrounded by a stabilizing membrane. Gene-based medicines are active at very low concentrations if they are effectively delivered to the target cells. Therefore, payload of drug in the microbubble is feasible for gene-based medicines with microbubbles. A number of studies have investigated the use of ultrasound without microbubbles in gene therapy [13]. In vitro studies have shown that ultrasound of sufficient intensity affects gene expression, up-regulating expression of some genes and down-regulating expression of other genes [14].
Nanovectors, particles with nano-scale dimensions are used for the delivery of therapeutic or diagnostic agents either through encapsulation or physical attachment of the desired moiety to the nanoparticle. Typically these systems are also composed of liposomes and many other macromolecular conjugates and particulate drug carriers. Stealth liposomes, the pegylated lipsomes are the prototypical example of first-generation nanovectors, whereby a nano-based delivery system enhances the delivery of a cytotoxic agent for improved therapeutic outcome. Stealth liposomes usually decorated by recognition molecules are being widely studied [15]. Nanoliposomes are also being currently used as they have improved therapeutic index due to the enhancement of permeability and retention effect [16].

Multiscale approach
Predicting the properties and behaviour of complex biological systems by computer simulation from a fundamental, ab initio perspective has long been a vision of computational scientists. It has been recognized since several years that the key to achieving this goal is utilizing hierarchies of methods and time scales that connect macroto nano-systems, i.e. a multiscale approach. Multiscale modeling and simulation has emerged as a new research area which has already had a significant impact on many scientific and engineering disciplines. Particular progress has been made in how to relate the molecular-scale chemistry to mesoscopic and macroscopic material and biomaterial properties. Based on large-scale atomistic and molecular modeling methods, hierarchical multi-scale modeling is capable of providing a bottom-up description of chemically complex materials [17][18][19]. Current perspective is to focus on liposomal delivery systems. The setting up of the multiscale methodology follows a simple procedure: starting from QM calculations on medium-size systems allows fitting parameters for the next atomistic step involving several thousand atoms, which then serves as a basis for coarse grained approaches. Despite the similarity of the strategy for different applications, the necessary parameters included in the approach are problem and system dependent. For example, the aim of this study necessitates describing a large area of a QM potential energy surface (PES), not only the minima, which will lead to different energy terms and different parameters. It is thus generally necessary either to proceed from step one, i.e. fit force field parameters on an adequately selected set of QM models, or, at least, test the available force field on the desired set of models and adjust it if necessary, before proceeding to the next step. Atomistic classical MM/MD calculations allow treating 50-60 Å systems in a time scale of several ns. Their use becomes rapidly unfeasible with the growth of the system dimensions. Even standard united atom approaches, by which methyl and methylene atoms are represented with a single Lennard-Jones term, might not be sufficient. For these reasons, coarse grain models have been proposed in the past few years, where a reasonable computational cost is reached through a reduction of the number of interaction sites [20]. There are several ways to make coarse graining, depending on the number of atoms which will be included in a grain (or bead) and depending on the nature and shape of the studied molecules. This approach has, in any case, to be based on a solid ground, i.e. on accurate enough atomistic calculations for the studied systems.

Force fields
Among the existing multiscale methodologies applied to phospholipid membranes, the first step, i.e. the fitting of MM parameters, has mostly been performed utilizing very small QM models, essentially for estimating the fixed atomic charges, the other parameters being adjusted along the time for a good description of experimental information [21,22]. The reason is that phospholipids are large molecules with a large number of degrees of freedom which could not be studied completely with accurate QM methods up to recently. Several phosphatidylcholine lipids, exploring their structures, conformations and dynamics [23] have now been studied.

Potential energy functions
The total energy of a biological complex system in classical mechanics is generally calculated as a sum of terms: (i) harmonic functions of bond lengths, bond angles, improper torsion angles; (ii) sinusoidal functions of torsion angles; (iii) nonbonding Lennard Jones interactions between atoms; (iv) electrostatic Coulomb interactions. This simple form, generally reported as Class I, allows treating very large systems of 100,000 atoms. Class II force field energy functions, include higher order terms to treat the bond and valence angle terms and/or cross terms between them, which increase the accuracy to treat conformational energies especially at geometries significantly far from the minimum-energy values. Further improvements may be performed with the explicit treatment of electronic polarizability. All MD studies available for lipid membranes have been performed with the simplest Class I potential energy function. Force fields for lipids can treat explicitly all atoms, including the numerous H atoms on the alkyl chains or neglect them to some extent, including them into a methylene or methyl group, or so called "united atom". In the case of lipid membranes, both all-atom and united atom force field have been used [24]. There are several possibilities to treat polarisability explicitly in the force field calculations. The most used methods are (i) using induced multipoles due to static ones [25,26], (ii) fluctuating charge models through redistribution of charges to reach equivalent electro negativity on each atom [27], (iii) fluctuating Drude model, where an additional Drude particle is attached to each atom and oscillates [28]; (iv) summing interactions between fragments ab initio calculated (SIBFA) [29]. In all cases, the polarizability has to be treated as a dynamical variable in the MD simulation, using extended Lagrangian methods [30]. To add up to the knowledge, no force field including polarization effects have been proposed up to now for lipids, except for the alkyl chains [31]. The majority of studies on polarizable force fields have treated water, both in gas and liquid phases, solvation of ions, ion -surfactant interactions, solvated small proteins or DNA [24].
The current perspective begins with focusing on QM/MM simulations followed by liposomal technology for controlled drug delivery systems. The methodology developed incorporates in the MM approaches, as much as possible, the information provided by QM studies of the properties explored at the mesoscopic level. Moreover, the effects of electronic intra-and intermolecular interactions, in particular the effects of polarization, has been included in the most approximate level of theory. An interface will be allowed to build by using different types of polarizable and reactive force fields, such as AMOEBA, GROMACS and CHARMM35. Such existing force fields will be first tested against QM results, in particular comparing MD simulations with respect to BOMD, and, then, adequately parametrized or reparametrized for phospholipids. Indeed, the evolution of a structure along time and with temperature has also to be properly mapped onto BOMD simulations. The dynamical behaviour of a single molecule is driven by the curvature of its potential energy surface, i.e. quantities similar to vibrational force constants, and by the vibrational energies. Intermolecular forces can be also tested comparing MD and BOMD results. Comparison of the different methods for treating polarization will lead to the choice of the most appropriate for lipid membranes in water. For the largest and most complex systems, full QM/MM computations are still too demanding. The use of the full MM method developed using the error-controlled strategy will be explored.

Methodology
The developed methodology will be used to study phospholipid membranes and vesicles, more particularly to the fluidity of the bilayers and its variation according to the nature of the lipids, addition of molecules inside the bilayer and presence of metal ions (Na + , Ca 2+ ) in the solution. As we all know the membrane is fluid at physiological temperatures and allows cells to change shape due to physical constraints or changing cellular volumes. The fluidity of the bilayer depends on the nature, topology and concentration of the various constituents. Fluidity is a vital property for the membrane to function. To compensate its variation with temperature, organisms use more unsaturated phospholipids or cholesterol. Studies of permeability in artificial bilayers indicate that fluidity and permeability of the membrane to small molecules may be related. The portion of the hydrocarbon chains adjacent to the headgroups is the likely site of the bilayer permeability barrier, selecting for solutes on the basis of partition coefficients, and acting as the major barrier to diffusion. Modeling the bilayer will indicate if this region of phospholipid layers shows a change in fluidity and if a reduced fluidity of the membrane may be reflected in a reduced ability of solutes to permeate the membrane.
Charged ions, or other larger molecules pass across the bilayer through channel compounds or ion carriers or, as recently proposed, water channels, but their presence at the bilayer surface may modify the dynamical properties of the phospholipids and therefore the membrane fluidity. Many biologically important properties of phospholipid bilayers are dependent on the pH of the environment, nature, and concentration of counterions. Among these properties are the phase transition temperature, packing, headgroup hydration, bilayer fusion, and formation of bilayer vesicles. Modeling the bilayer properties in the presence of ions allow quantifying the ionic sites in the bilayer and their impact on its physical dynamical properties.
The combination of important advances in methodology and computer codes and of new insight into the biological processes study will be able to make this perspective a solid and consistent ensemble, on top of which further levels of approximation can be built. Moreover, the error controlled strategy applied for optimizing an empirical method for phospholipids is novel in this domain. Combining these strengths in an approach that builds membrane models, integrating adequate atomistic and electronic information will represent a huge advance towards describing real membrane systems on a solid basis. The setup of the multiscale methodology follows a simple procedure i.e., starting from QM calculations on medium-size systems allows fitting parameters for the next atomistic step involving several thousand atoms, which then serves as a basis for coarse grained approaches.

Configuration and parameterization
Building phospholipid bilayer: The lipid bilayer was generated using a pre-equilibrated lipid membrane followed by energetic optimization and then the equilibrium run was carried out. Water molecules were specified by TIP3P parameterization and all lipids and water bonds were constrained with the SHAKE algorithm. Force field parameters were converted to GROMACS format and validation was ensured by comparison of individual energetic contributions to ensure proper porting. Equations of motion were inserted by using GROMACS simulation package. The initial bilayer system configuration consisted of a 16Χ12.5Χ 9.8 nm 3 box containing 46 DMPC molecules interspersed with 14 cholesterol molecules and solvated above and below with 1370 water molecules (Figure 1). The molecular composition of the leaflets is almost symmetric. The bilayer was parallel with the z-co-ordinate axes originating from the centre and extending 3.4nm into the aqueous phase. Periodic boundary conditions were used in all directions. vanderWalls and Coulomb cut-offs were set to 1.5nm, and the particle mesh Ewald summation was used for electrostatic interactions with the default associated parameters. Simulations were performed in the isothermal-isobaric ensemble at a temperature of a 308K and a pressure of 1 atm maintained by a Berendsen thermostat (τ=0.1ps) and Parinello-Rahman barostat (τ=1ps), respectively. Pressure coupling was applied isotropically to match the parameterization conditions. The sensitivity of the barostat scaling approach was tested by semi isotropic expansion in the xy-plane from that along the z-coordinate axis specifically designed for interface systems.
Potential mean force calculations: To calculate the free energy values along the bilayer and to avoid local hysteresis, the system has to be minimized to eliminate bad contacts (Figure 2 & 3). Further the system was equilibrated and calculations are being performed using constant force approach as implemented in the GROMACS pull code.   QM/MM deterministic approach has been applied where the built model tends to captures the short term dynamics and local structure. Also, at mesocopic level the model deals with the interaction of rough surfaces and long term dynamics with longer time scales. The thermal average over insertions of solute with randomly chosen orientations into configurations of the system at depth z unperturbed by the solute, the excess chemical potential, μ (z), and thereby the free energy of transfer, ΔG (z), from bulk water to the bilayer interior at the depth z were obtained by the following equation: ∆G (z) =μ (z) -μ (water) = -RTln<e -E(z)/RT >/<e -E(water)/RT >

Discussion
MD simulation studies suggest changes in the large scale properties of the bilayer but these changes are of local nature. However, the average effect on lipid chain ordering and thermotropic behavior is small. The formation of ordered domains involves co-operative interactions in terms of complexes, which in turn form networks through hydrogen bonding pathways. The current state is therefore very encouraging as it provides deep insight into the nature of atomic scale phenomena with a level of detail missing in any experimental technique. As a consequence, there are now fascinating issues that needs to be resolved such as the issues of domains in the multicomponent bilayers. The related effects on multicomponent bilayer due to interdigitation phenomenon and others need to be explored. Moreover, since the treatment of electrostatics is particularly important in lipid membrane systems, considerable attention needs to be paid to this issue. Two descriptions of partial charges are being employed. In the first case, no charges were considered and in the second approach, charges derived from ab inito quantum mechanical calculations using Gaussian with the Hartree-Fock method were used. The present model supports the computational data mining and also enables designing pro-drugs suitable for remote loading and reaching optimization of the loading conditions. Henceforth, we aim to propose a rational approach based on physicochemical principles which are expected to accelerate the development of novel liposomal formulations. The need to construct a quantitative model to predict the efficiency of drug loading into liposomes and the use of proper conjugation chemistry is mandatory. Finally, it is worth noting that the change in fluidity of the lipid bilayers, which is generated by the factors mentioned above, is a very important process, used to detect membrane alteration in tumors, thrombosis, etc. Also, conventional drug delivery systems have shown low efficiency and a continuous search for more advanced drug delivery principles is therefore of great paramount importance. High-performance computing has opened up new approaches for understanding phospholipids membrane models which may serve as molecular models for infectious diseases, metabolic disease and cancer. Theoretical model of loading drug into liposomes takes into consideration the multifactorial nature of the process. Liposomes are one of the best drug delivery systems for low molecular weight drugs, imaging agents, peptides, proteins and nucleic acids. The entry of targeted liposomal drug delivery onto the market is an innovative research idea within academia, government R&D that needs to be exploited in collaboration with the pharmaceutical industry. Drug Delivery Systems research is poised to become the tool to develop safe, stable, easy-to-use and, more importantly, economically affordable treatments.