Spontaneous Membrane Nanodomain Formation in the Absence or Presence of the Neurotransmitter Serotonin

Detailed knowledge on the formation of biomembrane domains, their structure, composition, and physical characteristics is scarce. Despite its frequently discussed importance in signaling, e.g., in obtaining localized non-homogeneous receptor compositions in the plasma membrane, the nanometer size as well as the dynamic and transient nature of domains impede their experimental characterization. In turn, atomistic molecular dynamics (MD) simulations combine both, high spatial and high temporal resolution. Here, using microsecond atomistic MD simulations, we characterize the spontaneous and unbiased formation of nano-domains in a plasma membrane model containing phosphatidylcholine (POPC), palmitoyl-sphingomyelin (PSM), and cholesterol (Chol) in the presence or absence of the neurotransmitter serotonin at different temperatures. In the ternary mixture, highly ordered and highly disordered domains of similar composition coexist at 303 K. The distinction of domains by lipid acyl chain order gets lost at lower temperatures of 298 and 294 K, suggesting a phase transition at ambient temperature. By comparison of domain ordering and composition, we demonstrate how the domain-specific binding of the neurotransmitter serotonin results in a modified domain lipid composition and a substantial downward shift of the phase transition temperature. Our simulations thus suggest a novel mode of action of neurotransmitters possibly of importance in neuronal signal transmission.


INTRODUCTION
Domain formation in biomembranes is thought to be essential for many biological processes as it may promote the localization of e.g., receptors in a specific sector of the membrane. Thereby, membrane domains facilitate as well the assembly or co-localization of different biological macromolecules by delimiting the accessible two-dimensional space on the membrane surface. Examples include immunological synapse formation (Grakoui et al., 1999;Davis and Dustin, 2004), or GPCR signaling related to dimerization Böckmann, 2016, 2020). Membrane domains are defined here as membrane regions confined in space and time, distinguished by membrane composition and/or physicochemical characteristics from their surroundings (Cebecauer et al., 2018;Kirsch and Böckmann, 2019). These domains have sizes ranging from a few nanometers covering in the order of 10-20 molecules only to more macroscopic domain sizes of a few hundred nanometers (see e.g., Cebecauer et al., 2018). In particular the direct visualization and characterization of small domains is experimentally limited to super-resolution methods such as stimulated emission depletion fluorescence correlation spectroscopy (STED) (Eggeling et al., 2009;Sezgin et al., 2017), interferometric scattering detection (iSCAT) (De Wit et al., 2015;Taylor et al., 2019), or atomic force microscopy (AFM) (Tokumasu et al., 2003;Shan and Wang, 2015). Further experimental evidence for small nanodomains was collected applying Förster resonance energy transfer (FRET) to model membranes (De Almeida et al., 2005;Petruzielo et al., 2013;Koukalová et al., 2017). Other nanoscopic information on membrane domains was achieved using X-ray scattering (SAXS, WAXS) (Sun et al., 2017), neutron scattering (SANS, QENS) (Pencer et al., 2007), or 2 H NMR (Veatch et al., 2007;Bartels et al., 2008;Bunge et al., 2008;Engberg et al., 2016;Bosse et al., 2019). However, while nanometer spatial resolution could be collected for a number of different model membrane compositions (see Cebecauer et al., 2018 for a comprehensive review), the sub-microsecond and typically as well sub-millisecond dynamics is masked by the limited temporal resolution of the above advanced microscopy techniques, or relies on model assumptions.
The gap in temporal resolution was in parts virtually closedand sometimes also blurred-by the progress in molecular dynamics (MD) simulations during the past two decades. In particular the so-called coarse-grained (CG) simulations characterized by condensation of the degrees of freedom by introduction of super-atoms or super-beads enabled the study of spontaneous domain formation in biomembrane models. The popular MARTINI model provided a first molecular view on domain formation of model membranes with fully saturated and polyunsaturated phospholipids as well as cholesterol, with well distinguished ordered ("raft" domain, or L O ) and disordered domains (L D ) (Risselada and Marrink, 2008) and has since then been applied in studies on domain formation for various scenarios (Friess et al., 2018;Lin et al., 2018;Bandara et al., 2019). Thereby, frequently little heterogeneity in the composition of individual domains emerged as a characteristics of coarse-grained simulations that is related to the coarsened lipid structure. That is, differences in lipids such as acyl chain length and degree of saturation may at CG resolution easily be exaggerated consequently resulting in pronounced domain separations. In addition, while CG simulations enable for long simulations (micro-to millisecond timescale) coupled to an enhanced dynamics of the membrane constituents (Marrink et al., 2007), the coarsening also affects membrane thermodynamics: The change in degrees of freedom results in a shifted balancing of enthalpy and entropy and thus possibly in a change of temperature-dependency e.g., in domain formation. Nonetheless, CG simulations are in many cases the method of choice to develop a molecular view on domain formation and domain structure.
The experimentally best characterized membrane domains are the sphingomyelin (SM) and sterol-enriched domains that compartmentalize cellular processes (Simons and Ikonen, 1997;Pike, 2006). The size of these domains initially coined membrane rafts was given as 10-200 nm (Pike, 2006). More recent experimental studies on model membranes containing phosphatidylcholine (PC), SM, and cholesterol at varying concentrations pointed to domain sizes between approximately 5 and 40 nm (Feigenson and Buboltz, 2001;London, 2011, 2015;Koukalová et al., 2017;Saitov et al., 2020). Frequently, these domains were discussed in the context of ordered (L O ) domains introduced early for phosphatidylcholine-cholesterol systems (Hjort Ipsen et al., 1987) and taken as model systems for rafts in cellular membranes. However, based on a combination of Monte Carlo simulations with FRET and z-scan FCS experiments, Koukalová et al. (2017) could show that SM and cholesterol driven nanodomains adopt a fluid and disordered L Dlike state with large amounts of PC lipids maintaining domain fluidity. These more dynamic domains were suggested as models of nanometer-sized heterogeneities or nanocompartments in biomembranes. A first atomistic-resolution picture of a SM-PC-Chol lipid bilayer from atomistic MD simulations was gained in 2015: A differential interaction pattern of Chol with SM and PC lipids was reported using experiment-derived compositions for L O and L D domains (Sodt et al., 2015). The L O -domain was described by a locally hexagonal substructure. Recently, also a spontaneous L O :L D phase separation was reported for DPPC:DOPC:Chol mixtures (Gu et al., 2020), accompanied by a significant partitioning of Chol to both ordered and disordered domains. In this study, the authors employed the Slipids force field (Jämbeck and Lyubartsev, 2012) shown earlier to perform well also close to the phase transition of PC bilayers (Pluhackova et al., 2016). A spontaneous phase separation or domain formation for sphingomyelin-containing PC membranes was, to our knowledge, not yet reported from atomistic simulations. According to Gibbs (1878), the word phase is used here synonymously for domains differing in composition or physical state observed for systems in thermodynamical equilibrium.
Membrane phase or domain formation is not only subject to changes in temperature, pressure, external fields, and lipid composition, but as well to (local) exposure to xenobiotics or endogenous molecules. Of particular interest is the interaction of neurotransmitters with membranes. These are released to the synaptic cleft in response to a signal. Their concentration in synaptic vesicles is as high as 270 mM (Bruns et al., 2000), suggesting a concentration of similar order within the synaptic cleft following exocytosis. Binding of serotonin (5-HT) to lipids would result in formation of a 2-dimensional neurotransmitter (NT) reservoir subject to 2D NT diffusion. Membrane binding drastically enlargens the receptor-neurotransmitter encounter rate and probably facilitates NT entry to membrane-buried receptor ligand-binding sites (Postila and Róg, 2020). A low pH and a high calcium level would in turn impede serotonin association within presynaptic vesicles (Mokkila et al., 2017). In a companion paper to this manuscript we describe by means of 2 H NMR spectroscopy that addition of serotonin to a POPC/PSM/Chol mixture induces increased differences in the lipid acyl chain order between ordered and disordered membrane domains . However, detailed knowledge on the interaction with and preference of serotonin (and other neurotransmitters) for plasma membrane (models) displaying L O and L D phases is lacking, albeit forming an important cornerstone in the understanding of the role of membranes for neurotransmitter dynamics and possibly linked neurological diseases.
The interaction of the neurotransmitter serotonin with simple phosphocholine lipid bilayers has been described before by means of atomistic MD simulations: Peters et al. (2013) analyzed the interaction of serotonin with single-component 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and 1,2dioleoyl-sn-glycero-3-phosphocholine (DOPC) for both charged and neutral forms of serotonin (pK a = 9.97 for primary amino group) on the 100 ns timescale. A protonation-dependent binding orientation of serotonin was reported, with the charged amine interacting with the lipid phosphate moiety and the deprotonated serotonin adopting a reversed orientation with the primary amine pointing toward the membrane core. The cationic serotonin form was independently used to investigate the effect of lipid acyl chain unsaturation (Azouzi et al., 2017): Serotonin was shown to preferentially bind to lipids with unsaturation on both chains. Here, however, both the hydroxyl group as well as the charged amino group of serotonin were seen to interact with the interfacial lipid headgroup region, in disagreement with the orientations reported by Peters et al. (2013), independent of the particular protonation state. Another simulation study addressed the binding of different neurotransmitters to model membranes mimicking intra-and extracellular membrane leaflets. The protonated form of serotonin was found to bind to both membrane compositions, albeit at a higher degree to the anionic intracellular membrane mimic (Postila and Róg, 2020). A charged serotonin was as well used in a simulation-based activation study of the 5-HT 3A serotonin receptor (Guros et al., 2020). Whether serotonin is indeed charged when bound to lipid membranes has to our knowledge not been studied. Also, little is known about differences in interaction patterns of neurotransmitters with ordered or disordered membrane domains, or lipid rafts.
Here, employing atomistic molecular dynamics simulations on the 10 µs timescale, we study the spontaneous formation of ordered and disordered domains for POPC/PSM/Chol mixtures at a molar ratio of 4:4:2 as a model for the extracellular leaflet of plasma membranes in the absence or presence of the neurotransmitter serotonin. The lipid composition was chosen to be in the L D /L O coexistence region at similar temperatures as studied here, as observed using 2 H NMR (Bartels et al., 2008;Engberg et al., 2016), EPR (Ionova et al., 2012), and fluorescence measurements (De Almeida et al., 2003;Veatch and Keller, 2005). Concomitant zeta potential measurements on vesicles of the same composition suggest that serotonin stays deprotonated, i.e., uncharged within the membrane. Overall, we observed an enhanced binding of serotonin to the L D domain resulting in changes in L D composition, decreased phase transition temperature, and enhanced dissemination of this phase.

Molecular Dynamics Simulations
We conducted 4-10 µs long atomistic simulations of POPC:PSM:cholesterol bilayers at a molar ration of 4:4:2 (mixed L O/D system), 8:61:31 (L O system), and 69:23:8 (L D system) in the absence or presence of 10 mol% serotonin (5-HT). The simulations were performed at a temperature of 303 K (L O , L D , L O/D ) and at 294 and 298 K (L O/D ). The lipid compositions for the L O and L D systems were selected according to a FRET-based phase diagram for a brain SM:DOPC:cholesterol composition using the tie line ends (Enoki et al., 2018) and a corresponding simulation study (Sodt et al., 2015). Similar tie line end compositions for L D and L O phases were reported for the ternary POPC:PSM:cholesterol system (De Almeida et al., 2003;Ionova et al., 2012). A summary of the simulated systems is provided in Table 1.
5-HT has two pH-sensitive groups: the primary amino group (NH 2 ) and an aromatic hydroxyl moiety (OH) ( Figure 5A), with pKa values of 9.97 and 10.73, respectively (Peters et al., 2013). Hence, serotonin in aqueous environment at neutral pH is cationic. However, the low dielectric constant environment (ǫ r ≈ 2 in the membrane hydrophobic core; Böckmann et al., 2008) inside the membrane could favor the neutral form of 5-HT, as suggested by the measured zeta potential of vesicles in the presence of 5-HT (see below). Since the binding preference to lipid bilayers was shown before for both the charged and the neutral forms of serotonin (Peters et al., 2013), we here modeled only the neutral form of serotonin as the assumed prevalent form of membrane-bound serotonin.
Bilayer systems were initially set up using the CHARMM-GUI web service (Lee et al., 2016). The total number of lipids in each system was between 508 and 510, and the lipids were solvated with at least 38 TIP3P waters per lipid, a salt concentration of 0.15 M NaCl, and an initial box size of ∼ 11 × 11 × 9.8 nm 3 . Serotonin molecules were added to the water phase of the systems at random orientation, and at a minimum distance of 2 nm from the lipid headgroups using the gmx insert-molecules tool implemented in GROMACS 2019 (Abraham et al., 2015). All the simulations were performed with the GROMACS 2019 software suite (Abraham et al., 2015) together with the CHARMM36 force field (Klauda et al., 2010). Serotonin parameters (Supplementary File 1) were generated with the CHARMM general force field (CGenFF) v4.0 (MacKerell et al., 1998;Vanommeslaeghe et al., 2010 using the ParamChem web server v2.2.0 (CGenFF, 2020).
Each system was minimized, heated, and equilibrated and then integrated under NPT conditions using a Nosé-Hoover thermostat (Nosé, 1984;Hoover, 1985), and a Parrinello-Rahman barostat (Parrinello and Rahman, 1981) with a compressibility of 4.5 × 10 −5 bar −1 . The integration timestep was set to 2 fs, and bonds to hydrogens were constrained with the LINCS algorithm (Hess et al., 1997). Long-range electrostatics were  L O/D computed using particle-mesh Ewald (Darden et al., 1993) with a tolerance of 10 −6 , 4th-order spline interpolation, and a maximal mesh size of 0.12 nm; van der Waals interactions were shifted smoothly to zero between 1.0 and 1.2 nm.

Data Analysis
Unless noted otherwise, the last 3 µs of each trajectory were used for the analyses, with frames taken every 1 ns. The analyses were performed using custom scripts involving the Scipy (Virtanen et al., 2020), Numpy (Oliphant, 2015), scikitlearn (Pedregosa et al., 2011), and MDAnalysis (Michaud-Agrawal et al., 2011) python packages. The HMM analysis was implemented based on the hmmlearn python library. Visualization and snapshot representation was done with VMD v 1.9.4 (Humphrey et al., 1996).

Analysis Protocol for Lipids Order States and Domains
In the L O/D coexistence systems, the state of the lipids was determined using a hidden Markov model (HMM) analysis, similar to the approach of Park and Im (2019). We decided not to use the local lipid composition as emission signal (Sodt et al., 2014;Gu et al., 2020) but rather the lipid order: Lipid chains and cholesterol orientation time series, quantified through the director order parameter P2 (Yankova et al., 2012), were used as observables to define three hidden states: putative L O , L D , and an intermediate ordered state. Briefly, the director order parameter is defined as P2 = 1/2 3cos 2 α − 1 , where α is the angle between the membrane normal and the end-to-end vector of the hydrophobic part of cholesterol or the hydrophobic tail of phospholipids. For phospholipids, |P2| was calculated as the average of both tails. The parameters of the model are the probabilities of each hidden state to have a certain P2 order, and the probabilities for a hidden state to stay or to change to the other one in a one-time step. The primary assumption is that the probability of the observable emission states, for each lipid type can be decomposed into those from high and lower order states and approximated by a combination of two normal distributions Park and Im, 2019). This was possible for all the lipid types in the simulated systems. As in Park and Im (2019), we trained a simple discrete HMM consisting of nine emission states, assigned by partitioning the observable space into nine subspaces (0,. . . ,8). The lowest and highest emission states are assigned when P2 < µ L and P2 > µ O , respectively. The intermediate transition states For the low and high hidden ordered states, the initial emission probability matrix was set up making use of the integrals of the normal distributions over the subspace of the observables. For the intermediate ordered state, we considered integrals over a normal distribution with variance to ensure that the intermediate states initial emission probabilities are localized within 3σ between the upper, µ O , and the lower, µ L , bounds (Park and Im, 2019). Given these initial set of parameters, the Baum-Welch algorithm (Baum et al., 1970;Rabiner, 1989) was used to find the maximum likelihood estimate of the parameters given as input the time series of the emission states for each lipid type in a given leaflet. Once the HMM parameters are determined, the most likely order state sequence for each lipid was determined by the Viterbi algorithm (Viterbi, 1967).
For better comparison of the hidden ordered, disordered, and intermediate states between simulations with/without serotonin and between simulations at different temperatures (L O/D systems), the emission states and the initial transition and emission probabilities were based on the POPC/PSM/Chol 4:4:2 simulation at 303 K.
Time-weighted histograms were computed to analyze the dynamical properties of the assigned ordered/disordered hidden states to POPC and PSM lipids. This analysis was performed on the last 3 µs of the L O/D system at 303 K with and without serotonin.

Calculation of Deuterium Order Parameters
Deuterium order parameters were calculated at each position along the aliphatic chains of POPC and PSM according to where θ is the angle between the C-H bond and the membrane normal (taken to align with z, while the bilayer is in the xyplane). Angular brackets denote an average over all sampled configurations. The order parameters were calculated by first averaging over time separately for each lipid in the system and then calculating the average and the standard error of the mean over the different lipids.
The lipids are divided into three populations for the calculation based on the HMM analysis.

Density Probability Distributions
Density probability distributions of water, lipids headgroups glycerol (POPC) and sphingosine backbones (PSM), and serotonin centers of mass along the bilayer normal, z, with respect to the bilayer center were computed for the last microsecond using the gmx density GROMACS 2019 tool (Abraham et al., 2015). The profiles of serotonin were not symmetrized as unsymmetrized profiles provide a useful check on convergence.

Potential of Mean Force (PMF)
Potential of mean force profiles of serotonin's center of mass were calculated from the last µs of the trajectory according to F(z) = −k B Tlnp(z). p(z) is the probability of finding serotonin at position z, and z is the position along the bilayer normal with respect to the membrane center, T is the absolute temperature, and k B the Boltzmann constant. The probability density p(z) was estimated using a Gaussian mixture model (Glodek et al., 2013;Bochicchio et al., 2018).

Serotonin-Bilayer Absorption and Serotonin Orientation
Serotonin absorption and association with the bilayer was determined using its hydroxyl oxygen atom z-position along the bilayer normal. The bilayer was divided into a headgroup and a hydrophobic core region based on the density probability distributions of selected groups of the phospholipids and the sphingolipids along the bilayer normal, as described above. Specifically, the headgroup region was assigned as the region between the choline and the glycerols z-positional probability distribution, the hydrophobic core as the region corresponding to the hydrocarbon chains.
Serotonin orientation in the bilayer-associated states was analyzed in terms of its indole group orientation. Two angles, θ and α were used to define the indole orientation according to a reference coordinate system depicted in Figure 5. The zaxis is perpendicular to the indole rings system, bisecting both the benzene and the pyrrole rings. The x-axis bisects the plane of the benzene ring, and the y-axis is orthogonal to the xand z-axes. θ is the angle between the z-axis and the bilayer normal, and α defines the angle between the x-axis and the projection of the bilayer normal in the x-y plane. A θ value of 90 o corresponds to an indole orientation with its plane on average parallel to the bilayer normal. α values of ∼ 120 o degrees and ∼ 310 o correspond to an indole orientation in which the hydroxyl oxygen atom is oriented toward the water phase. The percentage of serotonin residing within the hydrophobic core or the headgroup region was determined by calculating the relative number of contacts of serotonin's center of mass and the phosphate groups or the lipid acyl chains, respectively.

Area Compressibility (K A )
The instantaneous area per lipid A l (t) for the bilayers was calculated as the area of the simulation cell A divided by the number of lipids per leaflet. This assumes that undulations are small so the difference in projected and local areas is negligible. The area compressibility K A is evaluated from: Where A is the average total area, δA 2 is the mean square fluctuation, k B the Boltzmann constant, and T the temperature. The estimation of the standard error was performed using the block average method (Allen and Tildesley, 2017), with blocks of 10 ns (L D ) to 150 ns (L O ), corresponding to twice the correlation time of the time series fluctuations, estimated from the normalized correlation function of the average total area fluctuations.

Membrane Bending Modulus (K C )
The bilayer bending modulus, K C , was determined using the method developed by Watson et al. (2012), which allows reliable estimates of K C to be extracted from simulations of modestly sized boxes. This method analyzes the spectral thermal fluctuations of the lipid director vector fieldn. The lipid director is a vector pointing from the lipid head to tail and serves as a mean to quantify the lipid orientation. The theoretical prediction for the power spectrum of the longitudinal component ofn reads: Deviations from Equation (3) are observed only over wavelengths shorter than 3 bilayer thicknesses (Watson et al., 2012). Previous simulations have shown that simulations with 288 lipids only are sufficient to determine K C (Levine et al., 2014). For an in-depth description of the method, the reader is referred to Watson et al. (2012) and Levine et al. (2014).

Lipid Sample Preparation and Zeta Potential Measurements
Lipids were mixed in organic solvents at the molar ratio POPC/PSM/Chol 4/4/2 and evaporated at 40 • C in a rotary evaporator. Afterwards the lipids were hydrated at 40 • C for 30 min using K 2 PO 4 buffer (20 mM K 2 PO 4 , 100 mM NaCl, 0.1 mM EGTA pH 7.4 in Milli-Q water), freeze thawed and finally extruded through two 100 nm polycarbonate filters at 40 • C to make unilamellar vesicles (LUV). Lipid samples were diluted from a LUV stock solution to a concentration of 0.3 mM for the zeta potential measurements. Concentration of 5-HT in the samples varied between 0 and 10 mM. Serotonin was added externally from a K 2 PO 4 buffer (12.5 mM 5-HT, 20 mM K 2 PO 4 , 100 mM NaCl, 0.1 mM EGTA pH 7.4 in Milli-Q water). The samples were measured at 25 • C after 10 min incubation after externally adding 5-HT. Zeta potential was measured using a Zeta potential analyzer (Brookhaven) with 10 runs and 10 cycles.

RESULTS
The spontaneous formation of biomembrane domains was addressed in atomistic MD simulations for a POPC/PSM/Chol mixture (4:4:2, L O/D system) in the absence or presence of serotonin at three different temperatures (294, 298, and 303 K, each 8-10 µs). In addition, for comparison of the domain characteristics and quantitative binding analysis of serotonin, MD simulations of previously in silico characterized L O and L D domain compositions (Sodt et al., 2015) were employed as well. Convergence of the different systems was evaluated by monitoring the (overall) area per lipid and the number of lipid-lipid contacts for the different species (Supplementary Figures 1-3). The bilayer area equilibrated within ≈ 1 µs with a small drift to smaller areas in particular for the L O/D system at low temperature (294 K). Similarly, we observe small drifts in the number of lipid-lipid contacts as a measure of the formation of differently composed domains. Here, we thus focus on the initial steps of domain formation in ternary mixtures containing spingomyelin in the presence or absence of the neurotransmitter serotonin, domain composition, and domain influence on the mechanical membrane properties. The simulation system size and simulation length (10 µs) does not allow to conclude on the formation of large domains beyond 10 nm. Also domain compositions and sizes analyzed below may likely further slowly adapt on the 100 µs to millisecond timescale.
Membrane domain formation was analyzed using a Hidden Markov Model (HMM) employing P2 order parameterdistinguished emission states (see section 2 for details). For an unambigous assignment and for quantitative comparison of systems with/without serotonin as well as for comparison of simulations at different temperatures, the emission states were defined based on the L O/D simulation at 303 K.

Serotonin Charge
As a measure of the charge of serotonin bound to phospholipid membranes we determined the zeta potential for a POPC/PSM/Chol (4:4:2) mixture in the absence or presence of varying concentrations of serotonin. The zeta potential (see Table 2) was hardly affected even for high serotonin concentrations and varied between −0.1 and −4.5 mV, suggesting that serotonin preferentially binds in an uncharged i.e., deprotonatated form to the membrane. Also for pure POPC membranes no effect of serotonin on the zeta potential could be seen (data not shown). Accordingly, serotonin was chosen uncharged in all simulations.

Spontaneous Nanodomain Formation in Sphingomyelin-Rich PC Membranes
Membrane domain formation was monitored for the L O/D simulation (POPC/PSM/Chol 4:4:2 mixture). The HMM analysis revealed spontaneously formed ordered, sphingomyelinenriched domains (referred to as L O domain, see Figure 1 for snapshots) that include between 48% of all lipids (303 K) and 60% (294 K, compare Table 3). At 303 K, the L O domains contained ≈ 43% PSM, 34% POPC, and 23% Chol while the composition of the more disordered L D domain changed hardly with respect to the overall 4:4:2 composition (see Table 3). Interestingly, the L I domain (including lipids in intermediate ordered state) is observed to be enriched in POPC (52%). That is, the sphingomyelin-enriched L O domains are surrounded by PC-enriched intermediate regions that transition into a more mixed POPC/PSM L D domain. That indeed coexisting domains are formed is clearly discernible from representative snapshots displayed in Figure 2 (see also corresponding movies in the Supplementary Material): In particular at 303 K, larger regions of disordered domains formed that display a more loose lipid packing. In contrast, the POPC/PSM/Chol system at 294 K is mostly found in its ordered state (i.e., yellow) with only small more disordered patches, however at a comparably high lipid packing density. These differences in the formation of loosely packed L D domains are as well reflected in the average area per lipid. At 294 K, the area per lipid is ≈ 53 Å 2 and increases to ≈ 56.5 Å 2 at 303 K (Supplementary Figure 1, area per lipid excluding cholesterol). The observed spontaneous formation of domains differing in composition and physical characteristics suggests that the domains approximate corresponding coexistent phases in thermodynamic equilibrium.
In the presence of serotonin, the phase characteristics between L D and L O differed more significantly. Even at 294 K, a larger disordered domain is formed, lipids in disordered state now account for 26% of all lipids (18% w/o serotonin). Also, the domain composition changed upon addition of serotonin, in particular a significant change in the PC content is observed: While the PC and PSM content in the absence of serotonin is similar for the L D domain, addition of serotonin led to a significant enrichment of PC within the disordered domain at all investigated temperatures (e.g., from 39 to 48% at 303 K, Table 3). These phase changes are coupled to preferential binding of 5-HT to L D -phase membrane domains (triangles in Figure 2, lower panel) as discussed in more detail below.
To further assess the dynamical domain properties we analyzed the lifetime of POPC and PSM lipids within the   ordered and disordered hidden states for the L O/D system at 303 K with and without serotonin. In the absence of serotonin, POPC lipid lifetimes showed two clearly distinct populations for ordered and disordered state lipids (Supplementary Figure 4A). The time-weighted distributions were fitted by a double exponential and yield lifetimes of 17 and 51 ns for lipids in the disordered state as compared to 32 and 118 ns for lipids in the ordered state (fits performed on initial 300 ns). Addition of serotonin resulted in a significant stabilization of lipids within the disordered state, the lifetimes for the L D lipids increased to 20 and 99 ns (Supplementary Figure 4B). Note, however, that the lifetimes are expected to change for longer simulation times with increasing domain sizes. This effect is clearly seen in the number of disordered state lipids with lifetimes beyond 300 ns that is substantially increased in the presence of serotonin with increasing L D membrane domain sizes. In contrast, the lifetimes of PSM lipids in disordered and ordered states were not significantly affected by addition of serotonin on the short timescale (Supplementary Figure 5). Similar to POPC, the number of PSM lipids that were stabilized within the disordered state increased. These results are in line with the observed preferential binding of serotonin molecules to the L D phase of membranes (see below).

Phase-Dependent Lipid Order
The lipid acyl chain order provides information on the degree of ordering in ordered and disordered membrane domains. The deuterium order parameters for L O and L D phases as determined from the HMM differed substantially for the POPC/PSM/Chol (4:4:2) lipid membranes at 303 K, the average POPC sn-1-chain order |S CH | being 0.347 ± 0.008 for the L O phase and 0.224 ± 0.014 for the L D phase (see Figures 3, 4). The difference in lipid order diminishes for decreasing temperatures, the L D sn-1-chain order increases to 0.294 ± 0.017 at 294 K while the L O sn-1-order stays approximately constant (0.367 ± 0.006). The comparably high order for lipids assigned a disordered state by the HMM at low temperature is related to the small L D domain size and hence increased packing density of the lipids (compare Figure 2). That is, lipids in more disordered state are embedded in an overall highly ordered environment. A comparison to 2 H NMR at 303 K shows a substantially increased POPC sn-1 order for the L O/D system in the simulations . This discrepancy is probably related to an overall shift of the transition temperature to larger values observed for different atomistic force fields (Pluhackova et al., 2016).
The order parameters within the differently ordered domains of the mixed system at 303 K thereby agree well with the palmitoyl chain order determined for pure L D and L O systems with compositions equaling those of a previous simulation study (Sodt et al., 2014). Related to the different composition of these 1-phase systems-with POPC:PSM:Chol composition in L O at 8:61:31 and L D at 69:23:8-the difference between the average deuterium order parameters for L O and L D of POPC is comparable to the L O/D system with |S CH | O = 0.376 ± 0.009 and |S CH | D = 0.255 ± 0.042.
Upon addition of serotonin, a strong decrease of the lipid acyl chain order is observed for the disordered domains (for both POPC and PSM) at all studied temperatures. In contrast, the order within the L O domains was hardly affected by serotonin (see Figure 4). While in the absence of serotonin, truly disordered domains are observed only above 298 K, addition of serotonin results in large disordered and loosely packed domains even at 294 K. That is, the serotonin-lipid interaction partially disrupts the lipid packing. The increased disorder of L D domains as well as their increased size also affect the lipid area: At 303 K, the average lipid area is increased by ≈ 5% with respect to the serotonin-free case (Supplementary Figure 1).
Our results thus suggest a significant serotonin-induced decrease in the phase transition temperature for the POPC/PSM/Chol mixture.

Serotonin-Lipid Interaction
At simulation start, serotonin was added to the solvent phase of each system at random positions at a lipid:serotonin ratio of approximately 10:1 (in total 51 serotonin molecules, see Table 1). Spontaneous binding of serotonin to lipids was distinguished into binding to either the hydrophobic core region or to the  (Figure 5). In turn, 21 and 17% of all serotonin molecules bound to the hydrophobic core region of membranes with predefined L D composition or the L O/D membrane, respectively. Also within the L O/D membrane, we observed preferential binding of serotonin to the disordered membrane domains (compare Figure 1 and Figure 2).
The different adsorption modalities of serotonin concerning the membrane phases are well described by the position probability distribution profiles of its center of mass compared to that of selected lipid and cholesterol segments, and water, as well as by the potential of mean force (PMF) profiles. The profiles are shown in Figure 6 and the selected lipids groups are reported therein. The profiles for the phospholipid atoms and water molecules resemble the features that are generally observed in PC-rich bilayers, i.e., broad phosphate and choline distribution profiles, and penetration of water molecules into the lipid headgroup region including cholesterol hydroxyl groups. All the profiles clearly show that serotonin is able to penetrate the bilayer, irrespective of the membrane phase. For the L O -system, the highest probability for 5-HT is within the lipid headgroup region, that is, between the choline group and the glycerol/amine backbones. More than 50% of the serotonin remains in the solvent phase, only ≈ 3% are able to reach the carbon chains. An interesting density drop at the interface between water and the headgroup region is observed. It corresponds to a small but significant free energy barrier in the PMF profile (Figure 6) at the position of the choline group of PSM. In the PC-rich L D system, serotonin penetrates on average deeper into the membrane. For the L O/D phase the profile is broadened. The same tendency is indicated by the PMF data as well.
The generally broad 5-HT profiles reflect the mobility of the solute in the membranes. An analysis of the orientation of serotonin bound to the membranes is reported in Figure 5. We characterize the orientational distribution of 5-HT by attaching a local coordinate frame (Figure 5) to its indole group. The bilayer is symmetric under rotations about the bilayer normal, so two angles suffice to characterize the indole distribution: θ , the angle between the z axis of the indole and the bilayer normal; and α, the angle that is assumed when the bilayer normal is projected into the molecular x-y plane (compare section 2). Figure 5 shows the probability density of the indoles along these two axes. When bound to the headgroup region, serotonin shows broad orientational distributions for all the considered systems (right panel). However, orientations with θ = 50 o or θ = 135 o are preferred. This can be related to a tendency to form cation-pi stacking interactions with the choline headgroups, experimentally observed for indoles binding to PC lipids (Gaede et al., 2005).
The 5-HT orientation becomes restricted only when it reaches the glycerol backbone and enters the bilayer hydrophobic core (Figure 5, left panel). The orientation is dominated by structures Red plots (left) correspond to serotonin molecules residing in the hydrophobic core of the bilayer, whereas blue plots (right) correspond to those located within the lipid headgroup region. The percentages of serotonin residing at the hydrophobic core of the bilayer or within the headgroups is also reported.
with θ near 90 o -indicating that the molecular plane of the indole is orthogonal to the plane of the bilayer. This is particularly the case for serotonin binding to the L O system, while a much broader profile is observed for the L D and the mixed L O/D systems (see Figure 5). The same preferred orientation for the neutral form of serotonin was already observed by Peters et al. (2013) in a previous MD simulation study employing POPC, DOPC, and DPPC membranes. These maxima occur at α values of 120 and 300 o . In the former orientation, the hydroxyl group is oriented toward the water phase, and the aliphatic amino group points to the center of the bilayer. In the latter orientation, the molecule is rotated by 180 o corresponding to serotonin binding to the lower membrane leaflet (with reference to the coordinate system in Figure 5A). That is, for 5-HT embedded within the L O domains, the carbon chain aligns to the lipid acyl chains.
These results, together with the huge PMF barrier of > 7 kcal/mol in the membrane's interior (Figure 6), suggest a very small permeability for this solute, despite high partitioning into the membrane.

Membranes Mechanical Properties
The bending modulus K C is arguably the single most important quantity in membrane biophysics (Watson et al., 2012) and determines the membrane's ability to resist bending in a number of biologically relevant processes such as membrane fusion (Liu et al., 2006), membrane trafficking of molecules (McMahon and Gallop, 2005), and endocytosis (Chernomordik et al., 1995). As such, larger values of K C correspond to greater rigidity of the bilayers. The area compressibility modulus, K A , quantifies the response of the membrane area to tension, which under physiological conditions may arise from various perturbations including the addition of small molecules to the membrane.
To our knowledge, no experimental data-area compressibility and bending modulus-are available for the lipid bilayer compositions studied in this work. Hence, a direct comparison of our results to experiment is difficult. However, for the L D system (see Table 4) the average area per lipid (APL) is in good agreement with a previous simulation for the same composition (Sodt et al., 2015), however well below the APL reported for pure POPC bilayers in experiment (Kučerka et al., 2011) and simulations (Pluhackova et al., 2016). This decreased APL arises due to the condensing effect in particular of cholesterol. In addition, the area compressibility modulus, K A , and the bending rigidity, K C , are in line with previous studies of pure PC bilayers at 303 K, reporting K A = 180 − 330 mN/m (Binder and Gawrisch, 2001), and K C = 9 − 13 × 10 −20 J (Binder and Gawrisch, 2001).
Extremely large K A values for mixtures of PSM-Chol (1,718 mN/m) with respect to those of POPC-PSM mixtures (781 mN/m) have been experimentally reported by Needham and Nunn (1990), and related to an additional effect of PSM to decrease area fluctuations, possibly related to the tendency of SM to form intermolecular hydrogen bonds with cholesterol (Henriksen et al., 2004). In line with this observation, we observed a similarly large area compressibility of K A ≈2,140 mN/m for the PSM-rich L O system (61% PSM, 8% Chol). The average area per lipid for this system is also in full agreement with that reported for the same composition by Sodt et al. (2015). Our values for the bending rigidity of the L D system of K C ≈ 13.4 × 10 −20 J are in good agreement with experimental values for PC systems (Nagle, 2017). As the  Serotonin absorption and binding to the membrane has a strong effect on all the calculated physical observables. The data show a 16%, and 12 − 18% decrease in the bending modulus K C for the L D system and L O/D systems upon serotonin binding, respectively, as reported in Table 4. Such an effect, together with an increase in the average area per lipid and a strong decrease in K A , indicate substantial membrane softening upon neurotransmitter binding. Little or no change is observed for the pure L O phase (Table 4), as expected from the observed low serotonin absorption.

DISCUSSION
Plasma membranes are rich in phosphatidylcholines, cholesterol, and sphingomyelin. In this study, we provide an atomistic view on membrane domain formation of a ternary POPC/PSM/Chol lipid bilayer as a model for plasma membranes. Employing atomistic MD simulations, we trace the initial steps of the formation of (nano-)domains differing in composition, structure, and ordering i.e., dynamics, and mechanics, starting from a random mixture of the membrane constituents. Additionally, we demonstrate that a high neurotransmitter concentration may heavily affect the equilibrium between disordered and ordered lipid membrane domains by preferential binding to L D -domains.
In the ternary mixture, domains characterized by a substantial difference in lipid acyl chain order are formed only above 298 K. Interestingly, and despite the excellent agreement of these different acyl chain orders with those of simulations of systems with predefined compositions for ordered and disordered domains (see also Sodt et al., 2015), changes in domain composition are less distinct on the 10 µs timescale. That is, the ratios of PSM or POPC do not significantly differ between the differently ordered domains. This finding that composition lags behind ordering is most likely partly related to the limited size of the domains formed on the accessible timescale. It is important to note that ordered lipid domains (L O ) were defined in literature as domains containing lipids with extended and ordered acyl chains similar to the gel phase or L β -phase, and at the same time display a high mobility more similar to the disordered L D -phase (Brown and London, 1998;Polozov and Gawrisch, 2006). This rather counterintuitive characteristics, high order combined with high mobility, was supported by pulsed field NMR spectroscopy revealing a lipid diffusion coefficient that differs by only a factor of approximately 3-10 between ordered and disordered phases (Filippov et al., 2004;Orädd et al., 2005;Scheidt et al., 2005). Oppositely, experimental studies showed a comparable fluidity of ordered and gel phases using fluorescence probes (M'Baye et al., 2008). While the simulations yield a large acyl chain order, the lipid diffusion differs substantially between the ordered and disordered phases: The POPC diffusion coefficient within the L D -phase is 4.1 · 10 −8 cm 2 /s and is lowered by a factor of ≈ 25 within the ordered phase (diffusion coefficients analyzed for the predefined L O -and L D -phases, see also Pluhackova et al., 2016). Possibly, the different timescales between NMR (several ms) and MD (nanoseconds) impede the comparison. We speculate that, due to the long timescale, NMR may partly average the dynamics of lipids exchanging between ordered and disordered phases and thereby yield an efficient enhanced diffusion.
Most simulation studies dealing with serotonin assumed a protonated amine group of serotonin in agreement with the reported pK a of 9.97. However, binding of the neurotransmitter to lipid membranes with a low dielectric constant ( decreasing from ≈ 78 in the solvent to ǫ r ≈ 2 within the membrane hydrophobic core; Böckmann et al., 2008) as well as altered hydration properties may substantially shift the pK a (Narzi et al., 2008;Sandoval et al., 2016). A recent experimental study addressing the association of serotonin to model lipid bilayers revealed a nonspecific binding of the neurotransmitter at physiologically relevant concentrations (Josey et al., 2020). In the same work, neutron reflectometry results favored a model with deep penetration of serotonin into the bilayer, with its long axis oriented along the membrane normal. A protonated serotonin is at variance with the latter configuration since either the polar hydroxyl group or the charged amine would need to protrude into the hydrocarbon core region of the membrane. This indirect hint to a deprotonated membranebound form of serotonin is strongly supported by our zeta potential measurements on lipid vesicles displaying a largely serotonin-independent potential. Accordingly, our simulation study on the influence of serotonin on membrane domain formation focused on the deprotonated i.e., uncharged form of serotonin.
Addition of the neurotransmitter led to a significant downward shift of the phase transition temperature, thereby resulting in clearly distinguished ordered and disordered phases now even at 294 K. This phase shift is induced by preferential binding of serotonin to the domains containing mostly lipids in disordered state, enhancing thereby the fluidity of these domains. This is in agreement with a shift of the main phase transition temperature to lower temperatures reported before for DMPC upon addition of serotonin (Seeger et al., 2007) and similarly for ethanol (Griepernau et al., 2007). However, in the former study the shift induced by 3 mM serotonin was below 1 K. The L D domains under the influence of the bound neurotransmitters are characterized by an enhanced ratio of phosphatidylcholines, an increased area per lipid, as well as a drastic reduction in the bending modulus and the area compressibility. The latter findings are in apparent contradiction to a companion manuscript  reporting in particular 2 H NMR and AFM indentation experiments on the same ternary biomembrane models. In this work, an increased intendation force was required to rupture the membrane in presence of serotonin and interpreted as an increase in membrane stiffness. These opposing findings, a decreased bending modulus and area compressibility on the nano-scale and an increased membrane stiffness on the AFM scale, may be reconciled assuming that serotonin induces the formation of larger domains as supported by our 2 H NMR experiments . Additionally, serotonin induces a more pronounced disorder of the L D domains (this work and 2 H NMR; Engberg et al., 2020) as well as a more pronounced order of the L O domains ( 2 H NMR). The latter may be explained by changed domain compositions as observed in the simulations displaying an increased PC content of the L D domains and a thus likely decreased PC content for (larger) ordered domains and thereby increased ordering. Interestingly, in good agreement with our simulation data, a recent work reports a significant softening of lipid bilayers containing POPC/POPG/Chol (1:1:1) by serotonin .
In summary, the initial steps of the spontaneous lateral segregation of a plasma membrane model containing POPC/PSM/Chol (4:4:2) is observed on the microsecond timescale. Despite small differences in compositions, formed nano-domains display a huge difference in order resembling those for PC-rich disordered and PSM/Chol-enriched raftlike ordered domains. The neurotransmitter serotonin reshapes in particular the disordered membrane domains resulting in a PC-enrichment and a decreased membrane elasticity and probably in increased membrane domain sizes as suggested by comparison to 2 H NMR . These changes in membrane structure, domain formation and domain size, and membrane dynamics will probably among others affect receptor sorting and oligomerization in plasma membranes (Gahbauer and Böckmann, 2020), as well as membrane stability, and related permeability and pore formation Böckmann, 2016, 2019). A possibly similar mechanism was recently revealed for inhalation anesthetics that were shown to increase the size of ganglioside (GM1) rafts and to change the localization of the phospholipase D2 to GM1-enriched domains (Pavel et al., 2020).
Regarding its action as neurotransmitter, our results suggest a fast binding of serotonin to the postsynaptic membrane following synaptic vesicle fusion and thereby a fast and efficient mechanism for neurotransmitter removal from the synaptic cleft (Postila et al., 2016). The preferential binding of serotonin to L D -domains with a concomitant change in composition, ordering, and reduced phase transition temperature likely as well alters the membrane lateral pressure profile as reported before also for ethanol (Griepernau and Böckmann, 2008). A changed membrane pressure may in turn be hypothesized to shift the equilibrium between activated/non-activated receptor configurations if assuming different cross-sectional areas for different receptor activation states (Cantor, 1997a(Cantor, ,b, 1998. In addition, high concentrations of membrane-bound neurotransmitters would result in efficient binding to membrane-buried receptor binding sites.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
RB designed the research. ABo setup and carried out the simulations. ABo and ABr analyzed the simulations. OE performed the zeta potential measurements and analysis. RB wrote the paper with contributions from all coauthors. All authors contributed to the article and approved the submitted version.