Notes on the Treatment of Charged Particles for Studying Cyclotide/Membrane Interactions with Dissipative Particle Dynamics

Different charge treatment approaches are examined for cyclotide-induced plasma membrane disruption by lipid extraction studied with dissipative particle dynamics. A pure Coulomb approach with truncated forces tuned to avoid individual strong ion pairing still reveals hidden statistical pairing effects that may lead to artificial membrane stabilization or distortion of cyclotide activity depending on the cyclotide’s charge state. While qualitative behavior is not affected in an apparent manner, more sensitive quantitative evaluations can be systematically biased. The findings suggest a charge smearing of point charges by an adequate charge distribution. For large mesoscopic simulation boxes, approximations for the Ewald sum to account for mirror charges due to periodic boundary conditions are of negligible influence.


Introduction
Biomolecular membrane processes often take place on a microsecond scale, involving large interacting molecular ensembles with millions of atoms. Consequently, an atomic resolution modeling approach requires billions of integration time steps, with each time step accounting for all mutual atomic interactions based on approximate molecular mechanics force fields. Corresponding simulation runs need substantial computational resources, with simulation times often in the order of weeks or months [1]. In contrast, coarse-grained mesoscopic simulation techniques such as dissipative particle dynamics (DPD) considerably reduce the necessary number of interacting particles and allow much longer integration time steps on the picosecond scale, as soft particle-particle interactions replace their hard atomistic equivalents. As a result, mesoscopic simulations are orders of magnitude faster, with simulation runs of molecular ensembles representing millions of atoms for microseconds being completed within hours or days on standard multicore workstations [2,3]. Conversely, mesoscopic simulations imply a much lower resolution above the atomic level with only isotropic particle-particle interactions that average nonbonded interactions at the atomic scale-limitations that may prevent an adequate description of the molecular processes in question. Moreover, for an adequate representation of biomolecular membrane compounds or peptides and proteins, additional charged particles should be introduced, with long-range electrostatic interactions that superimpose the concurrent short-range interaction framework. However, hard electrostatic interactions between charged particles are alien to a soft interacting system and may lead to possible unphysical artifacts such as artificial ion pairing between differently charged particles [4].
This communication discusses the established approaches to integrating charged particles into a soft DPD context involving the interaction of cyclic peptides (cyclotides) with biological plasma membranes. Cyclotides exhibit collaborative membrane-disrupting activities on the microsecond scale [5]. Membrane lipid extraction as a specific mode of action [6][7][8][9] can be successfully investigated by DPD simulations [10], which, in particular, allow quantitative assessment of the dynamics of membrane disruption to characterize a specific cyclotide/membrane system for comparative purposes [11]. The agreement of quantitative lipid extraction within a cyclotide/membrane "sandwich" model (see details below) with experimental trends was demonstrated for more than two dozen cyclotide/membrane systems. Moreover, the linear additivity of cyclotide activity for cyclotide mixtures could be successfully modeled. In this work, the cyclotide/membrane "sandwich" model is utilized for an evaluation of the different approaches to charge treatment to promote mesoscopic modeling of biomolecular membrane processes.

Methods
Dissipative particle dynamics (DPD) is a mesoscopic simulation technique for isothermal complex fluids and soft matter systems. It satisfies Galilean invariance and isotropy, conserves mass and momentum, and achieves a rigorous sampling of the canonical NVT ensemble due to soft particle pair potentials that diminish molecular entanglements or caging effects. DPD is expected to show correct hydrodynamic behavior and to obey the Navier-Stokes equations [12][13][14][15][16]. DPD particle trajectories are guided by Newton's equation of motion mass and spatial position of particle i; t, time; F i , total force on particle i exerted by other particles j; F DPD ij , conservative soft repulsive DPD force on particle i exerted by particle j; F D ij , dissipative (frictional) force; F R ij , random force; F B ij , conservative harmonic bond force; F E ij , conservative electrostatic force. where the total force on a particle exerted by other particles consists of a conservative, a dissipative (frictional), and a random part. The opposing dissipative and random forces depend on each other and act as a thermostat conserving the total momentum and introducing Brownian motion into the system. The conservative forces comprise soft DPD particle repulsions that represent averaged nonbonding interactions of uncharged DPD particles as well as possible harmonic springs between bonded or electrostatic interactions between charged particles. DPD particles, in general, may be arbitrarily defined as fluid packets [13]. Molecular fragment DPD is a bottom-up variant [16][17][18][19][20][21][22] that identifies each particle with a distinct small molecule of molar mass in the order of 100 Da. Larger molecules are composed of adequate smaller molecular fragment particles that are bonded by harmonic springs to mimic covalent connectivity and spatial 3D conformations.
For charge treatment within the DPD framework [4,[23][24][25][26], the electrostatic potential between two charged particles can be neatly arranged as a product of a classical Coulomb term (denoted with index C), a charge distribution term (index D), and a splitting term (index S): U E r ij = t C r ij t D r ij t S r ij U E , electrostatic potential energy in reduced DPD units; r ij , distance between particles i and j in reduced DPD length units.
The Coulomb term is given by The dimensionless electrostatics coupling constant [4] evaluates to R cuto f f T K e, elementary charge; k B , Boltzmann constant; T, temperature; ε 0 , vacuum permittivity; ε r , relative permittivity; R cuto f f , DPD cut-off length. The temperature dependence of the relative permittivity for water can be approximated by [27] The charge distribution term turns the particle point charges into "smeared" charge distributions with a defined decay length [25] The splitting term (SP3 splitting in [26]) approximates the effect of mirror charges due to periodic boundary conditions (otherwise considered by the Ewald sum [23,24,26]): R cuto f f ,el , electrostatics cut-off radius in reduced DPD units.
The electrostatic forces on the charged particles are calculated by the derivatives of the electrostatic potential: ij F E ij , Electrostatic force on particle i exerted by particle j in reduced DPD units; r 0 ij , unit vector that points from particle j to particle i. If charge distribution and splitting are not taken into account, their corresponding terms are set to one (the pure Coulomb term remains, which will be denoted C); otherwise, the combinations of the Coulomb term with charge distribution (denoted CD) and splitting in addition (denoted CDS) are considered. Thus, the forces of the three different approaches evaluate to F E(C) ij for approach CDS.
For C alone, the resulting electrostatic forces are truncated to a maximum value to attenuate the hard potential by preventing its striving to infinity (where the maximum value is set to 25 reduced DPD units, similar to the repulsion of equal particles at room temperature with an electrostatic coupling constant of 1 as evaluated in [10] and used in [11]). All electrostatic interactions are confined to an electrostatic cut-off radius of 5 or 10 reduced DPD length units.
For studying membrane disruption by lipid extraction with molecular fragment DPD (with a single water molecule being represented by a single DPD particle), the versatile "sandwich" interaction model described in [11] is used where two plasma membranes surround an enclosed cyclotide/water compartment (compare Figure 3). The model itself is an artificial construct to estimate lipid extraction using a rapid simulation technique. Due to the lack of experimental values, the choice of cyclotide number in the cyclotide/water compartment is determined by maximizing the disruptive effect at the minimum cyclotide number, where the extremes would be a single cyclotide within the compartment (with possible membrane interaction but no disruptive effect) and a biologically unrealistic compartment consisting only of cyclotides without water. The model comprises the complete partitioning of the phospholipid molecules into molecular fragment particles, as well as a particle-based spatial 3D construction of the peptides with all DPD particle-particle repulsions where molecular particle-particle connectivity is controlled by additional harmonic springs. The rates of membrane disruption are determined as outlined in [11] with a minor difference: the evaluated rates describe the percentage of outer leaflet phospholipid molecules per microsecond that was extracted from the surrounding membranes as a more evident quantity in comparison to the percentage of ethane (Et) particles in [11]. Four cyclotide/membrane systems are analyzed that span the range of detected membrane disruption activities (see details in [11]): kB1-W19Y-P20S-V21T-L27T-P28S-V29T/NoC-PM with vanishing activity over kB1-W19Y-P20S-V21T/NoC-PM and kB1/NoC-PM to cO2-E6Q/NoC-PM with the highest activity. The acronyms denote the Möbius cyclotide Kalata B1 (kB1) and the Bracelet cyclotide Cycloviolacin O2 (cO2) with possible mutations, e.g., W19Y in kB1-W19Y-P20S-V21T denotes an exchange of hydrophobic amino acid tryptophan (W) in position 19 of wild-type Kalata B1 with the more polar amino acid tyrosine (Y). NoC-PM defines a biological plasma membrane that consists of a phospholipid composition of 40% DMPC, 20% DOPE, 5% PIP 2 , 10% DOPS, and 25% SM, with an inner to outer leaflet distribution of 24 to 76 for DMPC, 77.5 to 22.5 for DOPE, 75 to 25 for PIP 2 , 1 to 0 for DOPS (only inside), and 22 to 78 for SM, but without (uncharged) cholesterol. The chosen area per lipid leads to realistic membrane thicknesses and lipid distributions that are consistent with experimental findings and alternative simulation results. The simulation box size was doubled in the x-and y-direction (i.e., quadrupled in total) compared with the simulation boxes studied in [11], which led to an x-and y-dimension of 77 reduced DPD length units (corresponding to a physical length of about 377 Å) and a z-dimension of 104 reduced DPD length units (corresponding to a physical length of about 509 Å). Periodic boundary conditions are turned on in the x-and y-direction.
All simulations are carried out with the open DPD environment MFsim [28,29], which utilizes the open Jdpd simulation kernel [3,30]. The sketched electrostatic interactions are implemented in Jdpd classes ParticlePairElectrostaticsDpdPotentialCalculator (for the potential energy between two charged particles) and ParticlePairElectrostaticsDpdForceCon-servativeCalculator (for the electrostatic forces of charged particles) in method calculateParti-clePairInteraction, in which the methods could be easily extended to alternative calculation schemes (e.g., [31]). The MFsim graphical user interface is extended accordingly to allow for a comfortable setting of the electrostatic interaction parameters. All simulation job definitions are openly documented at [32] and may be viewed in full detail using the MFsim system. A simulation run of the cyclotide/membrane systems with nearly 2 million particles for 100,000 integration steps (corresponding to about 6 microseconds) performs within 15 h with 16 parallelized calculation threads. The calculation of the electrostatic interactions for charged particles requires less than 10% of the total simulation time with an electrostatic cut-off radius of 5 reduced DPD length units.

Results
The force functions of equally charged particles (with a valence of one) for the different electrostatic approaches are sketched in Figure 1. Since the splitting term depends on the electrostatic cut-off radius, the F CDS force functions are shown for the electrostatic cut-off radii of 5 and 10 reduced DPD length units, in which increasing the electrostatic cut-off radius leads to convergence of F CDS and F CD , with the splitting term approaching 1. While F C strives toward infinity for a vanishing particle-particle distance, F CD and F CDS run through a maximum toward a fixed force value that becomes zero for F CD .

Results
The force functions of equally charged particles (with a valence of one) for the different electrostatic approaches are sketched in Figure 1. Since the splitting term depends on the electrostatic cut-off radius, the FCDS force functions are shown for the electrostatic cutoff radii of 5 and 10 reduced DPD length units, in which increasing the electrostatic cutoff radius leads to convergence of FCDS and FCD, with the splitting term approaching 1. While FC strives toward infinity for a vanishing particle-particle distance, FCD and FCDS run through a maximum toward a fixed force value that becomes zero for FCD. Figure 1. Electrostatics forces F (in reduced DPD units) between equally charged particles (with a single elementary charge) for the different approaches (C, CD, CDS) along a distance relative to the DPD cut-off radius Rcutoff of one reduced DPD unit. Abbreviations are outlined in the text; the superscripts refer to the electrostatic cut-off radii of 5 or 10 reduced DPD length units.
In [10,11], a pure Coulomb approach was used to account for electrostatics particleparticle interactions with a truncated maximum force value and the electrostatic coupling constant being empirically evaluated to avoid artificial ion pairing. As a measure of the latter, the distance between differently charged ion pairs, which were initially located at the same position, was determined throughout the simulation (see Figure 1 in [10]). The evaluated electrostatic coupling constant led to distances of the differently charged ion pairs being equal to the corresponding uncharged particle pairs. While this approach may be plausible overall and prevents the initial ion pairs from simply sticking together, it neglects possible statistically strong ion pairing beyond enrichment or depletion of the charged particle coordination shells. To analyze the statistical surroundings of charged (valence 1) particles, a simulation run of 2 million particles (denoted as water H2O particles) with 50,000 positively charged (H2OP) and negatively charged (H2ON) ion pairs (H2OP-H2ON) was performed for 40,000 integration steps (corresponding to 2 microseconds). Figure 2 shows the resulting H2ON-H2OP and H2ON-H2ON radial distribution Figure 1. Electrostatics forces F (in reduced DPD units) between equally charged particles (with a single elementary charge) for the different approaches (C, CD, CDS) along a distance relative to the DPD cut-off radius R cutoff of one reduced DPD unit. Abbreviations are outlined in the text; the superscripts refer to the electrostatic cut-off radii of 5 or 10 reduced DPD length units.
In [10,11], a pure Coulomb approach was used to account for electrostatics particleparticle interactions with a truncated maximum force value and the electrostatic coupling constant being empirically evaluated to avoid artificial ion pairing. As a measure of the latter, the distance between differently charged ion pairs, which were initially located at the same position, was determined throughout the simulation (see Figure 1 in [10]). The evaluated electrostatic coupling constant led to distances of the differently charged ion pairs being equal to the corresponding uncharged particle pairs. While this approach may be plausible overall and prevents the initial ion pairs from simply sticking together, it neglects possible statistically strong ion pairing beyond enrichment or depletion of the charged particle coordination shells. To analyze the statistical surroundings of charged (valence 1) particles, a simulation run of 2 million particles (denoted as water H2O particles) with 50,000 positively charged (H2OP) and negatively charged (H2ON) ion pairs (H2OP-H2ON) was performed for 40,000 integration steps (corresponding to 2 microseconds). Figure 2 shows the resulting H2ON-H2OP and H2ON-H2ON radial distribution functions (RDF), which exhibit a distinct statistical ion pairing for the pure Coulomb approach in [10,11], whereas "charge smearing" by a charge distribution (F CD or F CDS ) leads to the expected statistical enrichment/depletion of the counter ion in the surrounding coordination shells. The RDFs for CD and CDS satisfy the relation g H2ON-H2OP (r) × g H2ON-H2ON (r) = g 2 (r), with g(r) being the RDF of uncharged H2O particles, which is consistent with the results of [4,25]. The different electrostatics cut-off radii of 5 and 10 reduced DPD length units lead to comparable results so that the shorter 5 reduced DPD length units (corresponding to a physical length of 22 Å) seem to be sufficient, which in turn leads to considerably decreased simulation times.
Membranes 2022, 12, x FOR PEER REVIEW 7 of 13 functions (RDF), which exhibit a distinct statistical ion pairing for the pure Coulomb approach in [10,11], whereas "charge smearing" by a charge distribution (FCD or FCDS) leads to the expected statistical enrichment/depletion of the counter ion in the surrounding coordination shells. The RDFs for CD and CDS satisfy the relation gH2ON-H2OP(r) × gH2ON-H2ON(r) = g 2 (r), with g(r) being the RDF of uncharged H2O particles, which is consistent with the results of [4,25]. The different electrostatics cut-off radii of 5 and 10 reduced DPD length units lead to comparable results so that the shorter 5 reduced DPD length units (corresponding to a physical length of 22 Å) seem to be sufficient, which in turn leads to considerably decreased simulation times. The cyclotide-induced membrane disruption by lipid extraction may be qualitatively characterized as follows [11]: From their initial random start distribution, the cyclotides begin to aggregate and form oligomers that eventually develop into tubular molecular superstructures. This collaborative cyclotide network allows membrane lipids to leave their membrane environment and distribute into the cyclotide/water compartment in between the plasma membranes (see Figure 3). In the course of time, the lipids of the outer membrane leaflets (directed toward the cyclotide/water compartment) are increasingly "consumed" by cyclotide superstructures, with lipid extraction becoming more and more "saturated." Membrane curvature or even rupture makes the quantitative evaluation procedure increasingly obsolete. Therefore, the evaluation procedure is limited to an intermediate "linear range" on the order of a few microseconds of the inherently nonlinear membrane disruption process. The cyclotide-induced membrane disruption by lipid extraction may be qualitatively characterized as follows [11]: From their initial random start distribution, the cyclotides begin to aggregate and form oligomers that eventually develop into tubular molecular superstructures. This collaborative cyclotide network allows membrane lipids to leave their membrane environment and distribute into the cyclotide/water compartment in between the plasma membranes (see Figure 3). In the course of time, the lipids of the outer membrane leaflets (directed toward the cyclotide/water compartment) are increasingly "consumed" by cyclotide superstructures, with lipid extraction becoming more and more "saturated." Membrane curvature or even rupture makes the quantitative evaluation procedure increasingly obsolete. Therefore, the evaluation procedure is limited to an intermediate "linear range" on the order of a few microseconds of the inherently nonlinear membrane disruption process. The details and the extent of the membrane disruption process are determined by the individual cyclotide and membrane types. The kB1-W19Y-P20S-V21T-L27T-P28S-V29T mutant, in which the important hydrophobic patch amino acids have been completely converted to more hydrophilic alternatives, still shows cyclotide oligomerization but does no longer form tubular superstructures that support membrane disruption. Thus, this mutant exhibits negligible membrane disruption activity. In contrast, the kB1 wild-type with an intact hydrophobic patch shows significant lipid extraction, while the kB1-W19Y-P20S-V21T mutant (half of the amino acids of the hydrophobic patch have been replaced by more hydrophilic ones) is in between (Figure 3). These qualitative findings are not affected by the different electrostatic approaches, i.e., the apparent order of activity is not changed.
However, the more sensitive quantitative evaluation of lipid extraction exhibits significant differences (see Figures 4-7). The pure electrostatic C approach leads to systematically reduced membrane disruption rates for all four cyclotide/membrane systems studied, whereas the CD and CDS approaches lead to comparable rates. A doubling of the electrostatic cut-off radius from 5 to 10 reduced DPD length units does not affect these findings. The pure C approach seems to overemphasize electrostatic interactions that lead to a specific membrane stabilization against disruptive attacks due to the charged phospholipid particles. This stabilization effect becomes more pronounced with increasing cyclotide activity for the kB1 variants, which all have the same zero net charge state (at pH 7.4). For the cO2-E6Q mutant, the stabilization effect is less pronounced in comparison to kB1, although its activity is higher. This may be traced to its different net charge state (+3 The details and the extent of the membrane disruption process are determined by the individual cyclotide and membrane types. The kB1-W19Y-P20S-V21T-L27T-P28S-V29T mutant, in which the important hydrophobic patch amino acids have been completely converted to more hydrophilic alternatives, still shows cyclotide oligomerization but does no longer form tubular superstructures that support membrane disruption. Thus, this mutant exhibits negligible membrane disruption activity. In contrast, the kB1 wild-type with an intact hydrophobic patch shows significant lipid extraction, while the kB1-W19Y-P20S-V21T mutant (half of the amino acids of the hydrophobic patch have been replaced by more hydrophilic ones) is in between (Figure 3). These qualitative findings are not affected by the different electrostatic approaches, i.e., the apparent order of activity is not changed.
However, the more sensitive quantitative evaluation of lipid extraction exhibits significant differences (see Figures 4-7). The pure electrostatic C approach leads to systematically reduced membrane disruption rates for all four cyclotide/membrane systems studied, whereas the CD and CDS approaches lead to comparable rates. A doubling of the electrostatic cut-off radius from 5 to 10 reduced DPD length units does not affect these findings. The pure C approach seems to overemphasize electrostatic interactions that lead to a specific membrane stabilization against disruptive attacks due to the charged phospholipid particles. This stabilization effect becomes more pronounced with increasing cyclotide activity for the kB1 variants, which all have the same zero net charge state (at pH 7.4). For the cO2-E6Q mutant, the stabilization effect is less pronounced in comparison to kB1, although its activity is higher. This may be traced to its different net charge state (+3 at pH 7.4), whereby the pure C approach may lead to an increased activity for this highly charged cyclotide that counteracts the membrane stabilization. at pH 7.4), whereby the pure C approach may lead to an increased activity for this highly charged cyclotide that counteracts the membrane stabilization.   Membranes 2022, 12, x FOR PEER REVIEW 9 of 13 at pH 7.4), whereby the pure C approach may lead to an increased activity for this highly charged cyclotide that counteracts the membrane stabilization.      Membranes 2022, 12, x FOR PEER REVIEW 10 of 13 Figure 6. Cyclotide/membrane system kB1/NoC-PM (see Figure 4 for details). Figure 7. Cyclotide/membrane system cO2-E6Q/NoC-PM (see Figure 4 for details). Figure 7. Cyclotide/membrane system cO2-E6Q/NoC-PM (see Figure 4 for details).

Discussion
Complex simulation models require numerous settings that follow theoretical considerations, deliberate choices, or simply experience. The influence of specific settings may vary for different areas of application. From a chemical point of view, the inclusion of electrostatic interactions is mandatory for an adequate description of biological membranes or biomolecules such as peptides or proteins.
In a DPD framework, a simple Coulomb approach with truncated forces tuned to avoid individual strong ion pairing still reveals hidden statistical pairing effects that may lead to artificial stabilization of molecular superstructures such as membranes or distortion of cyclotide activity depending on the cyclotide's charge state. While qualitative behavior is not affected in an apparent manner, more sensitive quantitative evaluations can be systematically biased. This was demonstrated in the complex disruptive interaction of cyclotides with a biological plasma membrane. The findings suggest a charge smearing of point charges by an adequate charge distribution. Since mesoscopic simulation boxes are large, there is only a negligible influence of mirror charges due to periodic boundary conditions. On the other hand, the comparable results of the CD and CDS approaches, in combination with their insensitivity to different electrostatics cut-off radii, show that it is advised to resist the temptation to overstretch the discussion of minor differences in model setups on the mesoscale. Adequate handling of charges may be particularly useful in studying the influence of a cyclotide's charge state on its membrane-disrupting activity, which is currently not well understood.  [29].