Potential energy function for a photo‐switchable lipid molecule

Photo‐switchable lipids are synthetic lipid molecules used in photo‐pharmacology to alter membrane lateral pressure and thus control opening and closing of mechanosensitive ion channels. The molecular picture of how photo‐switchable lipids interact with membranes or ion channels is poorly understood. To facilitate all‐atom simulations that could provide a molecular picture of membranes with photo‐switchable lipids, we derived force field parameters for atomistic computations of the azobenzene‐based fatty acid FAAzo‐4. We implemented a Phyton‐based algorithm to make the optimization of atomic partial charges more efficient. Overall, the parameters we derived give good description of the equilibrium structure, torsional properties, and non‐bonded interactions for the photo‐switchable lipid in its trans and cis intermediate states, and crystal lattice parameters for trans‐FAAzo‐4. These parameters can be extended to all‐atom descriptions of various photo‐switchable lipids that have an azobenzene moiety.

co-planar geometry of the two benzene rings in the crystal structure (Figures S1-S3) is consistent with earlier predictions from ab initio computation. [9] Understanding how azobenzene-based lipids interact with membranes, and how membrane properties change in the presence of azobenzene-based lipids, is important for the design of new photo-switchable lipids with tailored physical-chemical properties.
Presence of azobenzene-substituted pentyl phosphate (4-Azo-5P) in a membrane led to lower acyl lipid order parameters [5] ; the trans-isomer of 4-Azo-5P associated with larger anisotropy and compressibility of the area as compared to cis. [5] The trans-cis photoisomerization of photo-switchable lipids can be used to alter the local fluidity of a membrane [10] and the shape of unilamellar vesicles. [6] Computer simulations are a valuable approach to characterize the properties of lipid bilayers with photo-switchable lipids, because they allow us to model lipid membranes of various compositions and to explore the dynamics and molecular interactions of each of the components of the membrane system. The simulations involve solving numerically classical mechanical equations of motion for all atoms of the system, and then generating a molecular dynamics trajectory that describes the time-evolution of the atomic coordinates. Interactions between atoms of the system are described by a potential energy function that depends on the atomic coordinates of the system being studied, and on parameters-denoted force-field parameters, typically derived based on quantum mechanical (QM) computations. These parameters ensure that the structure and energetics of the system are properly described with the molecular mechanics (MM) force field computations.
MM force fields contain parameters for numerous types of lipid molecules. However, photo-switchable lipids are novel molecules that to date, lack an accurate atomistic force-field representation. Such a representation would require proper description of the azobenzene moiety, and of the bonded and non-bonded interactions between the azobenzene moiety and the remainder of the lipid molecule (Scheme 1).
The isolated azobenzene moiety and several azobenzene derivatives are represented in the Universal Force Field (UFF [11] ) and the Polymer Consistent Force Field (PCFF [12,13] ). These force fields can treat liquids and polymers, but are not optimized for simulations of biomolecular systems that include proteins. A coarse grain description of azobenzene was derived for MARTINI, [14] and a reactive force-field was derived for ReaxFF. [15] Several degrees of freedom of the azobenzene moiety were parametrized for AMBER: The dihedral angle that describes the twist of the benzene rings about the central N1 N2 double bond was parametrized using B3LYP/6-31G*, [16] and the C9 N1 N2 C3 and C9 N1 N2 dihedral angles (Scheme 1) were parametrized using Complete Active Space (CASSCF) computations. [17] A B3LYP-based AMBER parametrization for 4-hydroxy-4 0 -methyl-azobenzene [18] relied on atomic partial charges from electrostatic potential fitting (EPS). On the 4-OH-substituted benzene ring of this molecule, the charges for C and H atoms close to hydroxyl oxygen atom differ from those found on the unsubstituted benzene ring where C and H are close to 4 0 -methyl group. Description of azobenzene with alkyl-substitutions, as in FAAzo-4, might thus be difficult with this set of AMBER parameters.
Relatively recently, a set of CHARMM force field parameters was presented for azobenzene in McCullagh et al. [19,20] Atomic partial charges for azobenzene atoms were derived from ESP computations, [19] whereas for the internal consistency of CHARMM it is recommended that partial charges be derived by fitting to Hartree-Fock (HF) computations of water interaction energies. [21] Bonded parameters involving the nitrogen atoms (Scheme 1) were taken from parameters involving nucleic acid aromatic carbon, without further revision. [19] To enable reliable atomistic simulations of lipid membranes with photo-switchable lipids, we present here a complete set of CHARMM-compatible parameters for FAAzo-4 (Scheme 1). We used the CHARMM General Force Field (CGenFF) protocol [21] for the parametrization of drug-like molecules in CHARMM. In this protocol, MM charges are optimized by fitting to QM computations of water interaction energies and distances, and parameters for bonded degrees of freedom by fitting to QM computations of potential energy scans. [21] We obtained a good description of the geometry and water interactions of FAAzo-4 and of the crystal lattice from X-ray crystallography. The parameters we derived for FAAzo-4 can be transferred to the CHARMM parametrization of other azobenzene moieties with alkyl substitutions at positions 4-and 4 0 .

| Protocol used for deriving CHARMM force-field parameters
The empirical CHARMM MM force field equation [8] includes terms for bonded and non-bonded interactions as follows: where k b and k θ are force constants for the harmonic terms describing bond stretching and valence angle bending, respectively; b 0 and θ 0 are the equilibrium values for these terms. Dihedral angles φ are given as a sum of cosine functions with k φ -the force constant, n-the multiplicity, with integer values between 1 and 6, and δ-the phase shift. [8] Improper angles can be used to control chirality and planarity. [8] Nonbonded interactions between two atoms i and j separated by the distance r ij include electrostatic interactions between two atomic partial charges q i and q j , and van der Waals interactions characterized by the minimum interaction distance R min,ij , and the well depth ε ij. [8] In addition to the terms given by Equation (1), CHARMM includes the CMAP correction for the protein backbone. [22] To be consistent with the CHARMM force field, parameters for a non-standard molecule such as a photo-switchable lipid would need to be derived according to the CGenFF protocol for drug-like molecules. [21] This protocol involves QM geometry optimizations of the target molecule, optimizations of the atomic partial charges, and an iterative procedure to optimize bonded parameters. [21] The optimization of atomic partial charges relies on computations of water interaction energies and distances with QM and MM. Here, the water interaction energy ΔE is defined as the difference between the energy of an azobenzene molecule and a water molecule ( Figure 1c), and the sum of the energy E(azobenzene) of the isolated azobenzene molecule and that of the energy E(water) of an isolated water molecule: The water interaction distance R is the optimal interaction distance between the oxygen/hydrogen atom of the water molecule and the hydrogen or heavy atom of the azobenzene molecule ( Figure 1c).
ΔE and R are computed with internally fixed geometries of the azobenzene and water molecules. We used an azobenzene structure optimized with MP2/6-31G*; as recommended in the CGenFF protocol, [21] for the water geometry we used TIP3P. [23] The target QM values ΔE QM and R QM for all water-accessible donor and acceptor sites are computed first with HF/6-31G* [21] ; MM values ΔE MM , and R MM are then calculated, and partial atomic charges of the target molecule are adjusted manually until ΔE MM and R MM agree with the corresponding QM target values for all sites. [21] For neutral polar compounds, as is the case of the azobenzene moiety we parametrized here, ΔE MM values are scaled by 1.16, and R MM values are offset by −0.2 Å. [21] Optimization of the partial atomic charges is considered converged when ΔE and R computed with MM versus QM agree to within 0.2 kcal/mol and 0.2 Å, respectively. [21] For weak interactions, a less restrictive converge criterion of 0.5 kcal/mol can be used. [21] 2.2 | Identifying parameters of FAAzo-4 that require optimization To identify the CHARMM force field parameters that need to be refined for describing FAAzo-4, we used ParamChem [24] to search for parameters based on similarity to other molecules included in CGenFF. ParamChem searches for parameters and provides penalty scores that indicate the accuracy of the existing CHARMM parameters in describing the target molecule. A penalty score of zero indicates a good description of that particular degree of freedom, whereas a score >10 for a parameter indicates a poor force-field description.
Our ParamChem search indicated that the butyl and butanoic acid moieties of FAAzo-4 (Scheme 1) are well described by CGenFF, with penalty scores close to 0. Likewise, ParamChem indicated reasonable description for the partial atomic charges of carbon atoms distant from the central double bond (penalty score of 6 for C5, C7, C11, C13 in Figure 1a). In contrast, the partial atomic charges of atoms close to the central double bond (atoms C4, C8, C10, C14 in Figure 1a) had a score of 10, carbon atoms directly bonded to N1 and N2 a score of 127, and the partial atomic charges for N1 and N2 had a score of 139.
All bonded and non-bonded parameters involving N1 and N2 and the four carbon atoms directly bonded to the nitrogen atoms (C4, C8, C10, C14 in Figure 1a) had scores >10. A summary of the bonded parameters with high ParamChem scores is presented in Tables 1 and 2.
Pursuant to the observations above, we proceeded to optimize bonded parameters for the linkage region of 4,4 0 -dimethylazobenzene, and the partial atomic charges of the molecule. For the parametrization, we used the core structure of FAAzo-4, 4,4 0dimethylazobenzene (Scheme 1). For simplicity, in what follows we denote 4,4 0 -dimethylazobenzene as 4dmAB (Scheme 1, Figure 1a).

| Starting coordinates for 4dmAB in trans and cis isomeric states
The high-resolution crystal structure [4] was used for the starting coordinates of trans-FAAzo-4. The 4dmAB core moiety was prepared by replacing the sidechains of trans-FAAzo-4 with methyl groups using Avogadro. [25] The cis geometry of 4dmAB was prepared by editing z-matrix in Molden [26] to set the dihedral angles C9 N1 N2 C3 and N1 N2 C3 C8 ( Figure 1) to 8 and 56 respectively. [27]

| Geometry optimizations with QM
For all QM computations we used Gaussian09 [28] with standard convergence criteria. Our choice of MP2 for all QM geometry optimizations is based on tests indicating that MP2 allows us to reproduce the geometry of the cis-azobenzene fragment used for parametrization.
As a reference for the geometry of the cis-azobenzene fragment we used the crystal structure from Mostad and Rømming. [27] T A B L E 1 Summary of bonded parameters optimized for trans-4dmAB Atom names Note: The parameters are introduced in Equation (1). Atom labels are presented in Figure 1. All values reported in the table are MM-optimized parameters for bonded interactions.
T A B L E 2 Summary of bonded parameters optimized for cis-4dmAB Atom names 2.5 | Water interaction energy and distances and an automated protocol to optimize partial atomic charges We used the ForceField Toolkit [29] of Visual Molecular Dynamics (VMD [30] ) to generate starting geometries for the MP2-optimized structure of 4dmAB and water at 22 interaction sites ( Figure 1c). To calculate ΔE and R we used the parametrization scripts from the webpage of Prof. Alexander MacKerell Jr. at the University of Maryland, USA.
To make the charge optimization step more efficient, we implemented a Phyton script that uses CHARMM to compute optimal values for R and ΔE, and Sequential Least Squares Programming (SLSQP [31] ) in the SciPy Python library to optimize the partial atomic charges; the optimization of the charges was done in an outer Python loop that uses the multiprocessing Python library. The squared differences in ΔE and R obtained from the CHARMM calculation were used as a penalty function that is a subject to minimization with a convergence criterion of 1e −6 .
Test computations indicated that when 4dmAB is trans, atoms N1 and N2 ( Figure 1a) are poorly accessible to water, and we could not obtain convergence of the partial atomic charges for these two atoms.
For example, when we attempted to probe water interactions at N1 with the water molecule placed such that its hydrogen-bonding OH bond was perpendicular to the plane of the benzene ring, strong interactions between the water molecule and N2 hindered convergence of the QM computations. Likewise, when we probed water interaction such that the OH bond of the water molecule was parallel to the plane of the benzene rings ( Figure 1c) the repulsion from H7 or H12 atoms of 4dmAB would not allow the water molecule to hydrogen bond to N1 or N2 ( Figure 2, Table 3).
Since cis-4dmAB allows water to access the central double bond of the molecule and does not deviate strongly from trans-4dmAB charges (Table S1), we used this geometry to optimize all partial atomic charges.

| Choice of the van der Waals parameters
ParamChem assigned to the carbon and hydrogen atoms of azobenzene standard van der Waals parameters for benzene rings. In what follows, these standard van der Waals parameters are denoted as Set-1 ( Table 4).
All our attempts to optimize partial atomic charges of cis-4dmAB with van der Waals parameters Set-1 failed to give ΔE MM and R MM values within acceptable convergence criteria. To circumvent this problem, we used the van der Waals parameters used in CGenFF for aromatic ring hydrogen atoms adjacent to highly polarized heteroatoms such as nitrogen atoms. These van der Waals parameters are denoted here as Set-2 (Table 4).

| Optimization of the torsional potential for the central N1 N2 double bond
Test computations with MP2/6-31G* indicated that twisting the N1 N2 bond of trans-4dmAB to 50 costs 20 kcal/mol ( Figure S4).

| Derivation of parameters for bond stretching and valence angle bending
The initial bonded parameters obtained with ParamChem were tested and adjusted by performing three-point PES scans with MP2/6-31G*.
This procedure, which is consistent with the CHARMM force field protocol, [32] allows parametrization of the quadratic bond and angle potentials by fitting linear functions given by energies computed for three values of the degree of freedom being parametrized: The value of the degree of freedom (bond or angle) at the geometry-minimized structure without any restraint, one value decremented relative to this optimum, and one value incremented relative to the optimum.
Ideally, the incremented bond or angle values give energy differences of 1-3 kcal/mol relative to the energy-minimized structure. [32] PES profiles for bond stretching were computed with a step of 0.05 Å; for PES of valence angles we used a step of 4 . Bond lengths and equilibrium valence angles were fit to resemble MP2 or crystallographic structures.

| Crystal lattice simulation
To assess the usefulness of our parameters for describing equilibrium structures and inter-molecular interactions of FAAzo-4, we sought to use the parameters for simulations of the crystal unit obtained from X-ray experiments. [4] We first constructed coordinates for a crystal unit by using the crystallographic information file of the crystal structure of FAAzo-4 [4] and the free version of Mercury, the Cambridge Crystallographic Data Center crystal structure visualization software. [33] We then used these coordinates to generate a system of 64 trans-FAAzo-4 molecules using the Super Cell Builder functionality of Avogadro. [25] Crystal lattice simulations were performed using NAMD2 [34,35] using settings adopted from Nemkevich et al. [36] The time step for all simulations was 1 fs. We used a switching potential between 10 and 12 Å, particle mesh Ewald [37,38] for long-range Coulomb interactions, a Langevin piston to keep the pressure to 1 bar, and a Nose-Hoover thermostat to keep the temperature to the experimental value of 100 K; the Langevin damping coefficient was set to 1 ps −1 . Following geometry optimization and heating, we performed production runs of 2 ns using the NPT ensemble and isotropic pressure coupling.
In a separate set of simulations we constructed a crystal lattice composed of 64 4dmAB molecules based on the 4dmAB X-ray diffraction data from Harada et al. [39] The simulations on the 4dmAB crystal were performed using the same protocol as described above for FAAzo-4, but the temperature was set to 296 K, to be consistent with experimental conditions. [39] Convergence of the simulations was assessed by monitoring the crystal lattice vectors.

| MM potential of mean force profiles
To verify the dynamics of isolated cis-4dmAB at room temperature, we performed prolonged MM simulations from which we calculated potential of mean force (PMF) profiles for selected dihedral angles of the molecules. The geometry of the isolated cis-4dmAB molecule was first geometry optimized with the MM parameters derived here.
Molecular dynamics simulations were initiated by first heating the systems to 300 K, followed by equilibration with velocity rescaling for 500 ns, and then by a 5 μs production run. We used an integration step of 1 fs, and saved coordinates each 1 ps. From the simulation trajectory, we extracted the values of selected dihedral angles, and calculated histograms for these dihedral angles. We computed PMF profiles according to the equation where k B is the Boltzmann constant, T is the temperature in Kelvin, and density is the normalized distribution of the values of the dihedral angle.

| RESULTS AND DISCUSSION
We derived CHARMM force-field parameters for 4dmAB in trans and cis isomeric states and used these parameters to parametrize the com-

| QM geometry optimization of trans and cis-4dmAB
The MP2-optimized geometry of trans-4dmAB is very close to the starting crystal structure coordinates; the average root-mean-squared distance (rmsd) between the heavy atoms of the two structures is within 0.1 Å (Table 5). In the MP2-optimized structure, the N2 C3 bond has the same length as in the crystal structure, 1.4 Å; at 1.28 Å, the length of the N1 N2 bond in the MP2-optimized structure is within the 1.24-1.28 Å interval found in Böckmann et al. [41] from QM computations of azobenzene, and close to the 1.26-1.27 Å interval found from inspection of azobenzene-based molecular crystals. [39] The two benzene rings of the MP2-optimized trans-4dmAB structure are co-planar, which is compatible with previous QM computations on an isolated trans-azobenzene. [9] In the MP2-optimized structure of cis-4dmAB, the length of the N1 N2 covalent bond is overestimated by 0.2 Å relative to the crystal structure [39] (Table 6); the N2 C3 bond length is the same as in the crystal structure. [39] The C9 N1 N2 C3 dihedral angle is within 1 from the corresponding crystal structure value, the N2 C3 C4 C5 and N1 N2 C3 C4 dihedral angles within 3 ( Table 6).
The MP2-optimized geometry of cis-4dmAB is compatible with previous ab initio computations on a cis-azobenzene indicating a small twist of 0-11 around the central N1 N2 double bond. [9,41] 3.2 | Optimization of the partial atomic charges for transand cis-4dmAB Based on considerations of water accessibility to the central N1 N2 double bond of 4dmAB (Table 3), we used the cis isomer for all water interaction energy computations and optimization of partial atomic charges. We performed two separate optimizations by using different van der Waals parameters for hydrogen atoms of the benzene rings (see Set-1 vs. Set-2 in Tables 4 and S2). Our strategy was to derive partial atomic charges that fit the CGenFF convergence criteria for water interaction energies (within 0.2-0.5 kcal/mol) and water interaction distances (within 0.2 Å). [21] We initiated the optimization of partial atomic charges by including all atoms of the azobenzene moiety in the water interaction energy computations. As our test computations indicated poor convergence when attempting to re-compute all partial atomic charges, we reduced gradually the number of atoms whose partial atomic charges are allowed to vary; for atoms excluded from the charge optimization procedure we first used standard CHARMM partial atomic charges. These test computations led us to conclude that the CGenFF convergence criteria for water interaction energies are best met when only the partial atomic charges for N1, N2, C3, and C9 are optimized, and the partial atomic charges of atoms C6 and C12, which are zero in standard CHARMM, are set to 0.005e.
We found that the van der Waals parameters impact the description of the MM water interaction energies, and thus the convergence of the partial atomic charges. When using van der Waals parameters Set-1 (Tables 4 and S2) we obtained partial atomic charges that gave good description of ΔE and R for most of the heavy atoms of cis- Note: "MM initial" versus "MM optimized" indicate the starting CHARMM force field parameters versus the MM parameters optimized here. a MM simulations of a trans-FAAzo-4 crystal using the optimized MM parameters. b Crystal structure of trans-azobenzene. [40] c Crystal structure of trans-4dmAB. [39] d Crystal structure of trans-FAAzo-4. [4] e The value of 1.27 Å was estimated based on crystal structures of azobenzene compounds. [39] T A B L E 6 Geometry of cis-4dmAB optimized with MM versus MP2/6-31G*, and in the crystal structure  (Table 7).
To improve the description of the partial atomic charges of cis-4dmAB, we used van der Waals parameters Set-2 (Tables 4 and S2).
The agreement between ΔE QM and ΔE MM for all optimized sites improved to within −0.1 kcal/mol to 0.3 kcal/mol, and the R MM values computed for all interaction sites of cis-4dmAB were within the convergence criteria (Table 7). Moreover, the partial atomic charges optimized with Set-2 allow good description of the molecular dipole moment of cis-4DmAb ( Table 8): The dipole moment is overestimated by 17%, which is close to 20% overestimation recommended by the CGenFF protocol. [21] 3

.3 | Optimization of bonded parameters for trans-4dmAB
The central N1 N2 bond of the azobenzene ring was absent from the result of our ParamChem search for force-field parameters. As a reference bond length for N1 N2 of trans-4dmAB we used 1.27 Å, which is close to the 1.28 Å value in the MP2-optimized structure, and to the 1.26-1.27 Å values from crystal structures [39] ( Note: ΔE QM values were scaled by 1.16. [21] ΔE and ΔΔE values are reported in kcal/mol, and rounded up to the first decimal; and R and ΔR values in Å. Values in italics are larger than the CGenFF convergence criteria. Atom names are illustrated in Figure 1. Atomic partial charges optimized for Set-1 and Set-2 are presented in Table S3. computations a twist of 50 costs 20 kcal/mol ( Figure 4b). As we aim here to derive parameters to study the membranes with trans-or cis-FAAzo-4 close to the corresponding equilibrium structures, a description of highly twisted geometries close to the transition state for the isomerization is beyond the scope of our work. We thus restricted the parametrization of dihedral angles to torsions of up to 50 .
We obtained good agreement between the MP2 and MM PES for the C9 N1 N2 C3 dihedral angle by using two terms with multiplicities n = 2 and n = 4 ( Figure 4b, Table 5). We used the same multiplicities to describe the torsional profile of the N1 N2 C3 C4 dihedral angle that describes rotation of the benzene rings around the N C bond, and for the N2 C3 C4 C5 dihedral angle that describes bending of a benzene ring relative to the N1 N2 bond axis (Figures 4 and S4, Table 5). We found that the N2 C3 C4 H7 dihedral angle was well described once the torsional profiles mentioned above had been optimized, without the requirement of further optimization. The small tilt of 2 between the two benzene rings in the starting crystal structure [4] is likely due to inter-molecular interactions in the crystal. From MM simulations with the parameters derived here, we obtained an average tilt angle of 3.5 (Table 5).

| Force-field description of FAAzo-4
We used the parameters optimized here for 4dmAB to prepare a set of parameters for FAAzo-4 in its trans and cis isomeric states. As in our parametrization protocol we used standard CGenFF partial atomic charges for the benzene ring atoms and for methyl group atoms, we could replace the two methyl groups by butane and butanoic acid groups, respectively. Parameters for buthane and butanoic acid are transferred from CGenFF. The force field parameters for FAAzo-4 are presented in the Supporting Information.

| Crystal lattice simulation of trans-4dmAB and trans-FAAzo-4
To test the parameters we derived for trans-4dmAB, we performed crystal lattice simulations for trans-4dmAB at 296 K, and for trans-FAAzo-4 at 100 K ( Figures S1 and S2). The temperatures we chose are the same as used in crystallography.
Simulations of a crystal of 64 trans-FAAzo-4 molecules reproduced the size of the unit cell to within 0.3 Å with parameter Set-1, and to within 0.1 Å with Set-2 (Table 9, Figures S6 and S7).
Similarly small errors were obtained in simulations of the trans-4dmAB crystal at room temperature (Table 9).
We found that bond lengths, angles (Table 5), and inter-molecular distances ( Figures 5, S1c,d, and S2c,d) were accurately represented in the simulations. These computations, and the PES profiles discussed above, indicate that the parameters we derived here provide a reliable representation of trans-FAAzo-4.

| Optimization of bonded parameters for cis-4dmAB
In cis-4dmAB, two benzene rings are in close van der Waals contact ( Figure 6), and test computations indicated that we could not accurately describe structural properties of cis-4dmAB with parameters derived above for trans-4dmAB. To improve the description of the structural properties of cis-4dmAB, we increased the force constants for N1 N2 and N2 C3 bonds relative to trans-4dmAB, and decreased force constants for valence angles N1 N2 C3 and N2 C3 C8 (Table 2, Figures 6 and S8). In the one-dimensional PES for ζ, TS2 is sampled and its energy is the same, 1.7 kcal/mol, as that of TS2 in the two-dimensional PES. In the case of the ζ dihedral angle (Figure 9a), we could reproduce with MM the QM energy for TS1; for TS2, the MM and QM values agreed to within 1 kcal/mol (Figure 9a).
The MM-optimized structure of cis-4dmAB agrees well with the crystal structure of cis-azobenzene from Mostad and Rømming [27] ( Table 6) and with the MP2-optimized structure (RMSD = 0.07 Å between MP2 and MM). In the MP2-optimized structure and in the crystal structure, [27] the central N1 N2 bond is twisted by 8 . We obtained overall good agreement between MM with MP2 target values for bond stretching, angle bending ( Figures 6 and S8), and torsional profiles (Figures 8b, 9a, and S9).  (Figure 10a,b). At the MP2/6-31G* level of theory used for the parametrization, cis-4dmAB is 11.8 kcal/mol higher in energy than trans-4dmAB; with B3LYP, the energy difference is 15.5 kcal/mol (

| MM simulations and PMF computations of cis-4dmAB
The challenges we encountered when parametrizing dihedral angle terms for cis-4dmAB prompted us to use MM simulations to evaluate the dynamics of bond twists in the central region of cis-4dmAB.  To ensure the compatibility of the FAAzo-4 parameters with the CHARMM parameters for lipids and proteins, we used the CGenFF protocol for the parametrization of drug-like molecules. [21] The parameters we derived allow good description of the geometry and torsional properties of FAAzo-4. In molecular dynamics simulations of crystals, we obtain good agreement with X-ray diffraction data for the dimensions of the unit cells, molecular structure, and intermolecular distances (Tables 5, 9, Figures 5, S1 and S2).
The parametrization of the partial atomic charges was performed using cis-4dmAB, as for this conformer we could sample all sites in water interaction energy computations.
The partial atomic charges of cis-4dmAB allowed reliable param- The profiles for the water interaction energies of the central double bond for trans-4dmAB versus cis-4dmAB ( Figure 2) indicate that the cis photo-switchable lipid has stronger, that is, energetically more favorable water interactions than the trans lipid. Based on this result, we speculate that the trans photo-switchable lipid might locate deeper into the membrane plane than the cis lipid, whose double bond might rather prefer interactions closer to the lipid headgroup interface of the lipid bilayer.
The parameters derived here for the azobenzene moiety can be transferred to describe azobenzene-containing lipids for applications in photo-pharmacology-such as photo-switchable phosphatidylcholine, [6] diacylglycerol, [42] and ceramide. [10] As the parameters were derived using the CGenFF protocol for drug-like compounds, we suggest that they are applicable for simulations to probe, for example, the response of a lipid membrane to the presence of a small number of photoswitchable lipids, or the binding of a photo-switchable lipid to an ion channel. Atomistic models of lipid bilayers with high concentrations of photo-switchable lipids will likely require further developments from experiments and force-field parametrizations. Ana-Nicoleta Bondar https://orcid.org/0000-0003-2636-9773