Kai Steffen Stroh

Paulo C T Souza

Luca Monticelli

Herre Jelger Risselada
email: jelger.risselada@tu-dortmund.de

Paulo Cesar Telles De Souza

CGCompiler: Automated coarse-grained molecule parameterization via noise-resistant mixed-variable optimization

Keywords:

whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

are discrete variables whereas the bonded interactions are continuous variables. Here, we pioneer the utility of mixed-variable particle swarm optimization in automatically parameterizing molecules within the Martini 3 coarse-grained force-field by matching both structural (e.g., RDFs) as well as thermodynamic data (phase-transition temperatures). For sake of demonstration, we parameterize the linker of the lipid sphingomyelin.

The important advantage of our approach is that both bonded-and non-bonded interactions are simultaneously optimized while conserving the search efficiency of vector guided particle swarm optimization (PSO) methods over other metaheuristic search methods such as genetic algorithms. In addition, we explore noise-mitigation strategies in matching the phase transition temperatures of lipid membranes, where nucleation and concomitant hysteresis introduces a dominant noise term within the objective function. We propose that noise-resistant mixed-variable PSO methods can both improve as well as automate parameterization of molecules within building block CG FFs, such as Martini.

Introduction

Atomically detailed molecular dynamics (MD) simulations provide great insights into the structure and dynamics of biomolecular and other soft matter systems, but larger time-and length scales often require a coarse-grained (CG) description. In coarse-graining a group of atoms is mapped into one bead or supra-atom. Coarse-grained descriptions achieve computational efficiency by reducing degrees of freedom while preserving relevant aspects. This not only allows for bridging larger time and length scales but also enhances our understanding of the fundamental physics underlying molecular processes within biological cells. For example, it can enable fundamental insights into phenomena like the self-organization of lipid membranes and the formation of characteristic thermodynamic phases, including liquid-ordered, liquid-disordered, and gel phases. [START_REF] Risselada | The molecular face of lipid rafts in model membranes[END_REF][START_REF] Marrink | Simulation of gel phase formation and melting in lipid bilayers using a coarse grained model[END_REF][START_REF] Risselada | The freezing process of small lipid vesicles at molecular resolution[END_REF] Systematic coarse-graining approaches such as inverse Boltzmann and inverse Monte-Carlo approaches [START_REF] Lyubartsev | Calculation of effective interaction potentials from radial distribution functions: A reverse Monte Carlo approach[END_REF][START_REF] Lyubartsev | Multiscale modeling of lipids and lipid bilayers[END_REF] as well as force-matching approaches [START_REF] Izvekov | Effective force fields for condensed phase systems from ab initio molecular dynamics simulation: A new method for force-matching[END_REF][START_REF] Izvekov | A multiscale coarse-graining method for biomolecular systems[END_REF] parameterize coarse-grained force-fields by reproducing the structural part of the partition function of the fine-grained system by either matching relevant radial distribution functions or (combined) forces within the fine-grained system. However, because the partition function only describes a single thermodynamic state point at equilibrium, i.e., a unique combination of pressure & temperature values, systematically parameterized 'bottom-up' coarse-grained force-fields are not suited to describe phase transitions over a wider temperature range. Phase-transitions or phase-diagrams can, however, be optimally modeled using coarse-grained force-fields based on the alternative Statistical Associating Fluid Theory (SAFT) parameterization approach, which uses a scaled Lennard-Jones interaction potential whose functional form (the exponent) is uniquely adapted for each interaction type. [START_REF] Lafitte | Accurate statistical associating fluid theory for chain molecules formed from Mie segments[END_REF][START_REF] Papaioannou | Group contribution methodology based on the statistical associating fluid theory for heteronuclear molecules formed from Mie segments[END_REF] However, the main practical problem of all of these coarse-grained force-fields is their lack of chemical transferability, i.e. inclusion of a new molecule (interaction type) within the system would require reparameterization of all the existing interaction parameters. The Martini coarse-grained force-field [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF][START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] is a building block force-field, i.e., common chemical groups are parameterized as basic building blocks, which can be combined to build up any existing molecule. These basic building blocks of Martini, the beads, are parameterized top-down and reproduce the thermodynamic properties of the chemical groups they model, such as partitioning free energies in liquid-liquid systems, while complete molecules are parameterized with a combination of top-down (experimental data) and bottom-up (atomistic simulation). Such a parameterization enables the qualitative simulation of phase transitions as well as phase segregation in lipid membranes while simultaneously conserving molecular compatibility (transferability) by describing all non-bonded interactions with the same 12-6 Lennard-Jones potential form. However, a major drawback compared to other systematic coarse-grained approaches is that parameterization of molecules in Martini can be highly complex and often involves the selection and fine tuning of a large number of parameters (e.g., bead types and bond lengths) to optimally match multiple relevant targets simultaneously. A task that is time consuming when done by human labor. Additionally, it is not always obvious which parameters have to be changed in what manner to enhance a certain behavior, particularly when cooperative processes are involved. While the choice of individual bead types can be made using chemical intuition, still a sizable subset of combined possibilities exists. Importantly, parameterization of bonded and non-bonded parameters should be optimally performed simultaneously since bonded and non-bonded interactions are not independent -they are directly influencing each other via the density of interactions. [START_REF] Alessandri | Pitfalls of the Martini Model[END_REF][START_REF] Risselada | Martini 3: a coarse-grained force field with an eye for atomic detail[END_REF] Recent versions of the Martini force-fields such as Martini 3 rebalanced the density of interactions by introducing an even larger number of possible interaction types, thereby rendering the parameterization of molecules often a non-tractable problem to common users. Automation of coarse-graining is thus critical, especially when constructing large databases of molecules.Automation offers a solution to address the challenge of force-field development, which typically involves collaboration among multiple researchers working on interdependent parameters. By automating the process, a clear and structured flowchart-based hierarchy is established, providing an overview of how the parameterization is conducted and which objectives are targeted. This automation approach facilitates collaborations by allowing researchers to focus on selecting a set of relevant objectives and assigning importance or weights to each objective. These objectives, along with their individual weights, define the force-field's philosophy. Furthermore, automation empowers collaborations to prioritize two key aspects: the generation and provision of reference data for the objectives at hand, and the design of analysis tools to quantitatively assess how each objective is addressed within the automation pipeline. By automating the parameterization process, collaborators can allocate their efforts towards obtaining high-quality reference data that accurately represents the desired objectives. Simultaneously, they can focus on developing comprehensive analysis tools that enable thorough quantitative evaluation, ensuring the effectiveness of the automation pipeline in achieving the defined objectives. This collaborative approach maximizes the efficiency and reliability of the parameterization process while facilitating a deeper understanding of the force-field's performance.

Particle swarm optimization (PSO) is a powerful computational method used to optimize problems by iteratively improving candidate solutions based on a defined objective function.

Compared to evolutionary optimization methods like genetic algorithms, PSO offers advantages in efficiently finding global optima within high-dimensional continuous spaces due to its vectorial search direction. PSO has been successfully employed in various coarse-grained (CG) parameterization tasks, as demonstrated in previous studies. [START_REF] Bejagam | Machine-Learned Coarse-Grained Models[END_REF][START_REF] Bejagam | 4.0 opment of New Transferable Coarse-Grained Water Models[END_REF][START_REF] Empereur-Mot | Swarm-CG: Automatic Parametrization of Bonded Terms in MARTINI-Based Coarse-Grained Models of Simple to Complex Molecules via Fuzzy Self-Tuning Particle Swarm Optimization[END_REF][START_REF] Empereur-Mot | Automatic multi-objective optimization of coarse-grained lipid force fields using SwarmCG[END_REF][START_REF] Empereur-Mot | On the Automatic Optimization of Lipid Models in the Martini Force Field using SwarmCG[END_REF] PSO is primarily designed for continuous variables, making it well-suited for optimizing structure-based coarse-grained (CG) models where bonded and non-bonded parameters can be chosen from a continuum of values. However, in building block models like Martini, the non-bonded parameters are predefined and discrete, representing different interaction levels.

Consequently, the parameterization of molecules in a building block CG force field becomes a mixed-variable optimization problem.

When using PSO for parameterization in building block models, a transformation from the continuous space to the discrete space of force field parameters is necessary. This trans-formation introduces cumulative rounding errors, which can potentially affect the quality of the parameterization, especially in larger molecules. Therefore, additional evaluation and reparameterization steps are often required to ensure the optimal performance of the force field.

It is crucial to parameterize both bonded and non-bonded interactions simultaneously since they are not independent and their optimization should be performed in a coordinated manner. [START_REF] Risselada | Martini 3: a coarse-grained force field with an eye for atomic detail[END_REF] By considering their interplay during the parameterization process, the resulting force field can better capture the complex behavior of molecules in the system.

To address the limitations of existing PSO approaches, we employ a mixed-variable PSO scheme (mv-PSO) for parameterization. This approach allows for the simultaneous optimization of both discrete parameters (representing non-bonded interactions) and continuous parameters (representing bonded interactions), enhancing the accuracy and reliability of the parameterization process. Furthermore, due to the chaotic nature of MD simulations, observables measured in MD simulations are subject to noise. Since standard PSO was designed for deterministic objective functions, straightforward application to noisy optimization problems is error prone, because the algorithm can no longer correctly identify global and personal best solutions when noise levels are similar to differences between objective function values. 19 Noise-mitigation strategies are particularly important when utilizing thermodynamic data as targets, as these are notoriously expensive to estimate accurately in MD simulations, even when employing CG models. Particularly problematic is the targeting of phase transition temperatures, which involve a first order phase transition and are thus subject to nucleation and concomitant hysteresis.

In this paper, we pioneer the application of mixed-variable particle swarm optimization in automated parameterization of molecules within the Martini 3 coarse-grained force-field by matching both structural (e.g., RDFs) as well as thermodynamic data (phase-transition temperatures). The important advantage of this approach is that both bonded-and non-bonded interactions are simultaneously optimized while conserving the search efficiency of vector guided particle swarm methods over other metaheuristic search methods such as genetic algorithms. In addition, we explore noise-mitigation strategies in matching the phase transition temperatures, where nucleation and concomitant hysteresis introduces a dominant noise term within the objective function. To the best of our knowledge, the impact of noisy objective function values has not been previously addressed in the context of applying PSO for CG parameterization. The manuscript is structured in the following way: Section 2 describes the mixed-variable PSO algorithm and parameterization procedure. As an example, we parameterized the linker region of sphingolipids, a biological highly relevant class of lipid molecules, that constitutes approximately 30 mol% of the plasma membrane lipids, [START_REF] Van Meer | Membrane lipids: where they are and how they behave[END_REF] but has not been updated for Martini 3, yet. Details of the simulated molecules, systems and observables are given in Section 3. Results are presented in Section 4, followed by conclusions in Section 5.

ticle swarm optimization

With CGCompiler we present a Python package that streamlines CG molecule parameterization. It employs mixed-variable particle swarm optimization to simultaneously optimize categorical (beadtype) and continuous (bonds, angles, dihedrals, ...) variables. Therefore, CGCompiler is particularly well suited for, but not limited to, parameterization tasks in CG FFs that follow a building block approach. To enable the application of the building block approach also to larger molecular fragments, consisting of more than one CG bead, the method allows for optimization of shared building blocks in different molecules, e.g. the headgroup, linker, or tails of lipids.

Molecule parameterization in Martini 3 follows three steps: i) Choice of mapping and bead sizes ii) Assignment of chemical bead types iii) Choice of bonded terms and assignment of bonded parameters. [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] While a mapping from atomistic to CG model and the set of of bonded terms have to be predefined, the here-presented algorithm optimizes bead size, chemical bead type and bonded parameters simultaneously.

The parameterization workflow is shown in Figure 1. For a given parameterization task, the user provides or generates the target data, and creates a set of CG training systems, that allows measurement of the target observables. In the initial iteration, the optimization algorithm generates a number N p , i.e., the swarm size, of candidate solutions with random FF parameters, and runs MD simulations for each candidate solution and each training system.

Candidate solutions are then scored by how well the parameterization targets are reproduced.

By utilizing the swarm's knowledge of the fitness landscape, candidate solutions are updated, and a new cycle of MD simulations, analyses, and fitness evaluations, starts. This is repeated until a termination criterion is fulfilled. Due to noise in the objective function evaluation, the selection of the true best parameters can only be done with a certain probability. Therefore, the set of the best, statistically equal candidate solutions undergoes a screen-to-the-best procedure, which either provides one solution that is significantly better than the rest, or reduces the field of viable candidate solutions further, on which more expensive evaluation simulations would be performed.

Figure 1: Parameterization workflow. i) A set of training systems from which the target properties can be extracted. ii) Target data is acquired from atomistic simulations and experiments. iii) An initial swarm is generated with FF parameters randomly selected from a predefined range of feasible parameters. iv) All candidate solutions are simulated in all training systems, the target observables are measured and compared to the target data, i.e., the fitness of the candidate solutions is estimated. New candidate solutions are generate by utilizing the swarm's knowledge of the fitness landscape. v) Step iv is repeated until a termination criterion is fulfilled. vi) A screen-to-the best procedure yields the optimized set of FF parameters.

Mixed-variable particle swarm optimization

In the original PSO algorithm for continuous optimization problems in a D-dimensional parameter space, particle i has a position vector [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF] At each iteration t the velocity and position are updated by

X i = (x 1 i , ..., x D i ) and a velocity V i = (v 1 i , ..., v D i ).
V i (t + 1) = w * V i (t) + c 1 r 1 (pbest i (t) -X i (t)) + c 2 r 2 (gbest(t) -X i (t)) (1)
X i (t + 1) = X i (t) + V i (t + 1) (2)
Where pbest i (t) is the personal best position of particle i and gbest(t) is the best position found by the whole swarm. w is an inertia weight, which balances global vs. local search.

The coefficients c 1 and c 2 are balancing personal vs. social experience. r 1 and r 2 are vectors of random numbers. In the mv-PSO algorithm, that is utilized in our work, the position vector of a particle takes a hybrid form, where Z dimensions encode continuous variables and V dimensions encode categorical variables. 22

X i = (x 1 i , x 2 i , ..., x Z i continuous , x Z+1 i , x Z+2 i , ..., x Z+V i categorical ) (3)
The continuous and categorical parts of the position vector are updated separately.

Continuous reproduction method

In classical PSO the swarm can get trapped in local optima and therefore prematurely converge. [START_REF] Wang | A particle swarm optimization algorithm for mixedvariable optimization problems[END_REF] To promote diversity while maintaining good convergence efficiency Wang et al.

proposed an altered continuous reproduction scheme, where particle i learns from the best position of a randomly selected particle. [START_REF] Wang | A particle swarm optimization algorithm for mixedvariable optimization problems[END_REF] In order to guide the swarm towards improved solutions, the pool of pbest to choose from, only consists of solutions whose fitness is superior to pbest i (t).

V i (t + 1) = w • V i (t) + c • r • (pbest r (t) -X i (t)) (4)
Algorithm 1 Continuous reproduction method 1: Input: sorted swarm, particle i, parameter w i 2: for j = 1..Z do

3:

Randomly choose r, i ≤ r ≤ N 4:

v j i (t + 1) = w i • v j i (t) + c • r • (pbest j r -x j i )
5:

x j i (t + 1) = x j i (t) + v j i (t + 1) 6: end for 7: return (x 1 i , x 2 i , ..., x Z i )

Categorical reproduction method

Values of categorical variables are assigned according to a probability. Initial probabilities are given by

P rob j,n (0) = 1 n j (5)
where n j is the number of available values for the jth variable. To leverage the swarm's knowledge of good solutions, only the superior half of the sorted swarm is utilized in updating the probabilities of available categorical values. To avoid premature extinction of available values, a lower limit is assigned for P rob j,n . If P rob j,n falls below that lower limit, P rob j,n is set to that threshold value, and all probabilities are renormalized such that n P rob j,n = 1.

The categorical update method is shown in Algorithm 2.

Algorithm 2 Categorical reproduction method

1: Input: sorted swarm, particle i, parameter α i 2: for j = 1..V do

3:

for each available value n, n = 1 to n j do Count j,n = 0 4:

for each personal best pbest i , i = N/2 to N do 5:
if pbest i,j == V alues j,n then 6:

Count j,n + = 1 7: end if 8:
end for 9:

P rob j,n (t + 1) = α i • P rob j,n (t) + (1 -α i ) • Count j,n N/2 10:
end for 11: end for 12: for j = 1..V do

13:

Assign an available value to x Z+j i according P rob j 14: end for 15: return

(x Z+1 i , x Z+2 i , ..., x Z+V i )
Molecule parameterization is typically a multiobjective optimization problem (MOP). A simple way to scalarize an MOP is by linear weighting. The scalarized optimization problem is solved by minimizing the cost, which is given by

cost = o w o f o (x) (6)
Where w o is an objective weight, f o the objective cost function, and x the parameter vector.

The objective weights can be used to balance the importance of the utilized parameterization targets. The weights are set by the user. Setting weights might require some intuition about the parameterized molecule, quality of target data, etc.

Each objective can have a different objective cost function f o . New objective cost functions can be added by the user easily. In its present form, the parameterization algorithm uses two distinct objective cost functions. For single valued observables, such as area per lipid, membrane thickness, melting temperature, solvent accessible surface area (SASA) the objective cost function is defined as

f o (x) = 1 Ns s w o,s Ns s w o,s 1 N types,s Ntypes,s t max(0, SAE(y s,t (x), ŷs,t ) -E tol o,s ) . (7)
y s (x) is the observed value, given the FF parameters x. ŷs is the target value. N s is the number of training systems that is used for the current parameterization objective. N types is the number of bond or angle types being parameterized. The deviation from the target is calculated by the scaled absolute error SAE(y, ŷ) = ŷ-y ŷ . With the error tolerance E tol o,s , uncertainties in target data can be accounted for. Each training system has an additional weight w o,s , which can be used in case of differences in target data quality or similar cases.

Generally these are set to 1.

For observables that are given in the form of distributions, such as bond lengths, angles, or radial distribution functions (RDF), the objective cost function is given by:

f o (x) = 1 s w o,s s w o,s 1 N types,s Ntypes,s t EM D(φ(x s,t ), φs,t ) (8)
Where φ(x) is the observed distribution, given the FF parameters x. φ is the target distribution. The earth mover's distance EM D(φ(x s,t ), φs,t ) is a measure of the distance between the two distributions. [START_REF] Rubner | The Earth Mover's Distance as a Metric for Image Retrieval[END_REF]

Noise mitigation strategies for PSO

PSO was designed for deterministic objective functions. Due to the chaotic nature of MD simulations hereby measured observables are subject to noise. With noise in objective functions, selection of the true best solutions is not guaranteed. Since solutions, that are identified as the best, attract the swarm toward regions of interest in parameter space, noise can misguide the swarm and therefore deteriorate PSO performance.

Resampling

Resampling is a widely applied strategy for noise mitigation within the objective function.

Relatively simple resampling methods are equal resampling (PSO-ER), extended equal resampling PSO-EER, and equal resampling with allocation to top-N solutions PSO-ERN. [START_REF] Rada-Vilela | Population statistics for particle swarm optimization: Resampling methods in noisy optimization problems[END_REF] These simpler methods are regularly outperformed by state-of-the-art resampling methods, such as optimal computing budget allocation PSO-OCBA, [START_REF] Pan | Particle swarm optimization for function optimization in noisy environment[END_REF] but the quality of results depends on the specific optimization problem and noise levels. 19,[START_REF] Rada-Vilela | Population statistics for particle swarm optimization: Resampling methods in noisy optimization problems[END_REF] OCBA aims to maximize the probability of correctly selecting good solutions. This is done by first allocating a primary computational budget equally to all current solutions to estimate their cost means and variances. A secondary budget is then sequentially allocated to solutions with lower means and higher variances to improve the fitness estimations of potentially good solutions. For efficient secondary budget allocation at least 5 primary evaluations should be executed for mean and 13 https://doi.org/10.26434/chemrxiv-2023-r97bm ORCID: https://orcid.org/0000-0001-5239-7124 Content not peer-reviewed by ChemRxiv. License: CC BY-NC 4.0 variance estimation. [START_REF] Chen | Simulation Budget Allocation for Further Enhancing the Efficiency of Ordinal Optimization[END_REF] This might make application of OCBA prohibitively expensive for regular CG molecule parameterization tasks. Based on the observation that most observables utilized in the multiobjective optimization of the sphingomyelin linker region have a low variance and only a few suffer from a larger variance (cf. Figure S5), we hypothesize that in the molecule parameterization task at hand, one primary objective function evaluation is sufficient to differentiate potentially good solutions from bad solutions, but to maximize the probability of correctly selecting the true best solution, the accuracy of the fitness estimates has to be increased. Therefore, we propose a somewhat pragmatic approach, that salvages the core idea of OCBA, i.e., allocate additional computational budgets to where it is the most useful (low mean and high variance). At each iteration, our resampling method involves one full objective function evaluation of the current solutions. The current solutions are then ranked by their fitness, and for the best N solutions only the observables that have significant variance are reevaluated.

Set of statistically equivalent solutions

Even with noise mitigation, at the end of an optimization run, there will be a number of solutions with very similar scores. While in a deterministic setting, the global best position is determined by

gbest = arg min x∈Pt f (x), (9)
where P t is the set of all positions that have been visited by the swarm up to iteration t, with noise in the objective function no solution can be declared the best with 100% certainty. 19 With the screen-to-the-best procedure of Boesel et al. [START_REF] Boesel | Using Ranking and Selection to "Clean up" after Simulation Optimization[END_REF] a set of positions P g t ⊆ P t can be selected, such that the true global best solution gbest is contained in P g t with probability of at least 1 -α (with 0 < α < 1). 19 For solutions i, j ∈ P t , fi and S 2 i denote the sample mean and sample variance of objective function values. The elementary steps of the screen-to-the-best procedure are:

1. Compute W ij , W ij = t i S 2 i n i + t j S 2 j n j 1/2
, ∀ i = j ∈ P t (10) where

t i = t (1-α) 1/|P t |-1 ,n i -

Simulation details

The Python package is based on evo-MD. [START_REF] Methorst | Inverse design of cholesterol attracting transmembrane helices reveals a paradoxical role of hydrophobic length[END_REF] All simulations were performed with GROMACS 2020.4 and 2021.4 29 and analyzed with in-house Python scripts that are utilizing MDAnalysis, [START_REF] Michaud-Agrawal | A toolkit for the analysis of molecular dynamics simulations[END_REF][START_REF] Gowers | A Python Package for the Rapid Analysis of Molecular Dynamics Simulations[END_REF] LiPyphilic, [START_REF] Smith | A Python Toolkit for the Analysis of Lipid Membrane Simulations[END_REF] SciPy, 33 and pyemd, which is a Python wrapper for Pele and Werman's EMD implementation. [START_REF] Pele | Fast and robust Earth Mover's Distances[END_REF][START_REF] Pele | A Linear Time Histogram Metric for Improved SIFT Matching[END_REF] Visualization was done with NGLview. [START_REF] Nguyen | NGLview-interactive molecular graphics for Jupyter notebooks[END_REF]

Atomistic models

All atomistic models were simulated using the CHARMM36 37-39 force field. Table 1 provides details about the atomistic target systems. Initial configurations of the membrane systems were generated with the CHARMM-GUI membrane builder. [START_REF] Jo | CHARMM-GUI Membrane Builder for Mixed Bilayers and Its Application to Yeast Membranes[END_REF][START_REF] Wu | CHARMM-GUI Membrane Builder toward realistic biological membrane simulations[END_REF][START_REF] Lee | CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field[END_REF] Following energy minimization and equilibration, all systems were simulated with a 2 fs time step. Bonds of hydrogen atoms were constrained employing the LINCS algorithm. [START_REF] Hess | LINCS: A linear constraint solver for molecular simulations[END_REF] Van der Waals forces were gradually switched off between 1.0 nm and 1.2 nm. The PME algorithm [START_REF] Darden | Particle mesh Ewald: An N • log(N ) method for Ewald sums in large systems[END_REF]

Coarse-grained models

All coarse-grained models were simulated using the Martini 3 11 force field. Beta version 14 of the Martini 3 cholesterol parameters was used. [START_REF] Borges-Araújo | Martini 3 Coarse-Grained Force Field for cholesterol[END_REF]48 Initial configurations of membrane systems were generated with the Python script insane. [START_REF] Wassenaar | Computational Lipidomics with insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations[END_REF] Details of the employed training systems are listed in Table 2. All systems were energy minimized and equilibrated with the current version of DPSM, that made the Martini 2 model of sphingomyelin compatible with Martini 3. During the particle swarm optimization each system was equilibrated with the candidate FF parameters in two stages, with time steps of 2 fs and 20 fs, respectively. For all coarse-grained production simulations a time step of 20 fs was used. Non-bonded interactions were cut off at 1.1 nm. For electrostatic interactions the reaction-field method was used with a dielectric constant of 15 and the reaction-field dielectric constant was set to infinity.

Temperature coupling was obtained via the velocity rescale algorithm [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] with a coupling time τ t = 1.0 ps. System pressures were held at 1 bar by using the Parinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] with a coupling time τ p = 12.0 ps. Pressure coupling was applied isotropically for aqueous solutions and semi-isotropically for membrane systems. In simulations for melting temperature estimation anisotropic pressure coupling was employed, using the Berendsen barostat 50 with a coupling time τ p = 4.0 ps. Our aim was the development of an automatization framework for molecule parameterization in building-block force fields. As an example we parameterized the sphingolipid linker region.

Section 4.1 shows the results of the parameterization with CGCompiler using a simple noisemitigation strategy. Since noise-mitigation strategies can only reduce the effects of noise when selecting the true best solution, the best statistically equivalent solutions generated during the mv-PSO run are subsequently screened-to-the-best, as described in Section 2.2.2.

Parameterization of the sphingolipid linker region

Table 3 shows the observables and their weights used in the parameterization. The swarm size was 64. Noise-mitigation was done by reevaluating the melting temperature of the 16 best candidate solutions of the current iteration 12 times, i.e., results were obtained with noise-mitigation setting mv-PSO-R16 (cf. Section 4.2). As T m is the major contribution to cost variance, but the employed T m estimation method is good for differentiating good from bad solutions, i.e., it has an accuracy of a few K. Other observables were only evaluated once, area per lipid (APL) fluctuations were the second largest cause of cost variance. For more details on noise-mitigation efficacy see Section 4.2.

All results shown include the complete set of the best statistically equivalent candidate solutions P g that remained after two rounds of the screen-to-the-best procedure (cf. Section 2.2.2). This set contains 18 candidate solutions. to SM(18:0). [START_REF] Venable | CHARMM All-Atom Additive Force Field for Sphingomyelin: Elucidation of Hydrogen Bonding and of Positive Curvature[END_REF] It is therefore not unexpected that the Martini DPSM models show a reduced thickness compared to SM(18:0).

While the melting temperatures estimated with the biphasic approach, that is used during optimization for performance reasons, are not within the specified tolerance regime of 2 K but ≈ 5 -6 K below the target value and ≈ 3 -4 K below the lower target threshold, the new models are greatly improved compared to the current model, which was 20 K off target.

Notably, the estimation of T m is approach dependent. Estimations using the alternative, reversible melting approach with slow melting rates, based on Kowalik et al. [START_REF] Kowalik | Combination of MD Simulations with Two-State Kinetic Rate Modeling Elucidates the Chain Melting Transition of Phospholipid Bilayers for Different Hydration Levels[END_REF] and Sun and Böckmann 52 (see SI for further details), which requires a very large computational budget (as done here, total simulation time for one T m estimation > 90 µs) show an even better agreement with the experimental melting temperature.

The here-performed biphasic approach utilizes a bilayer that is half gel and half liquid. The gel phase is fabricated by quenching to a temperature well below the melting temperature, and the gel phase system is combined with a preequilibrated liquid system.

The combined system is then equilibrated with thermostats set to different temperatures for the two phases. As quenching and equilibration can take up to several hundreds of ns, reconstructing the starting structure for every candidate solution would significantly increase computational cost of a PSO run. Therefore, starting structures for this procedure were generated with the current DPSM parameters beforehand and equilibrated using the parameters of each candidate solution. While equilibration of the fluid phase is generally fast, this certainly is not the case for the gel phase. Considering that an unequilibrated phase is inherently less stable, the presence of an equilibrated liquid phase alongside an unequilibrated gel phase may lead to a slight systematic underestimation of the melting temperature (Tm). [START_REF] Coppock | Determination of Phase Transition Temperatures for Atomistic Models of Lipids from Temperature-Dependent Stripe Domain Growth Kinetics[END_REF] However, this potential underestimation can be anticipated and taken into account during the analysis.

The equilibrium melting rate approach does not suffer from this potential problem of unequally equilibrated phases. To minimize bias caused by the quenched starting structures used in this approach, for each validated candidate solution eight different starting conformations were generated. The solvent accessible surface area (SASA) is commonly used to further compare the molecular volumes and shapes between CG and AA models. [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF][START_REF] Borges-Araújo | Improved Parameterization of Phosphatidylinositide Lipid Headgroups for the Martini 3 Coarse-Grain Force Field[END_REF] Figure 5 shows the SASA values of P g in comparison to the AA and current CG DPSM models. The SASAs are computed for the linker beads AM1 and AM2, as well as all supra-atoms that are directly connected to the linker, i.e., beads PO4, T1A, and C1B, as these connections are also parameterized. With SASA values of ≈ 6.24 nm 2 all newly parameterized CG models show a better reproduction of the AA value (5.24 nm 2 ) compared to the current model (6.45 nm 2 ), but with discrepancy of ≈ 19 % all SASA values remain grossly too high. It appears that solely reparameterizing the linker region is not enough to fix this issue. Furthermore, using SASA directly as a target in the high-throughput optimization scheme is not necessarily beneficial, since a specific SASA value is not a unique representation of a certain shape.

Therefore, comparisons of solvent accessible surface areas between AA and CG models are most helpful when done by simultaneous visual inspection. For automated parameterization, however, more detailed shape descriptors should be used.

Force field parameters

Non-bonded interactions: Due to the polar nature of the linker region of sphingolipids, only the chemical types of the P-block of the Martini 3 FF were eligible. As groups of 3 or 4 heavy atoms were combined into supra-atoms in the specified mapping, bead sizes small (S) and regular (default) could be chosen by the algorithm. Both bead sizes were permitted for both interaction sites, to allow for some wiggle room, even though 4 heavy atoms are grouped together into supra-atom AM1 and 3 into AM2. A slight miscount of mapped atoms is not uncommon in Martini, e.g., the mapping of the NC3 bead is actually 6-to-1. [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF] One feature of the mixed-variable approach is that the optimization procedure directly yields a probability distribution of bead types, cf. Figure 6A. While for the interaction site AM2 there is clear consensus on the bead type, for AM1 only the size (small) is clearly determined, but there is some ambiguity regarding the interaction strength. The reduced size of one of the beads seems to be warranted, given the still too high SASA values shown above, and is also inline with the new Martini 3 models of glycerolipids. [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] It is also worth mentioning that the chemical bead types chosen by our algorithm match the expected assignment suggested by Martini 3.

A converged "degenerate" probability distribution of bead types is the result of two or more bead types having indistinguishable effects on fitness. This can be caused by noise levels being larger than the fitness differences or the employed set of observables and training systems is lacking the necessary discriminatory power. Both issues can be remedied in post-optimization screening, but should optimally be addressed during optimization. As the former option would merely improve selection the pool of generated candidate solutions, the later would potentially allow the generation of truly better solutions.

Additionally, for both, non-bonded and bonded FF parameters, diversity can be caused by the fact that the objective cost function for single valued observables (Eq. 7) has an error tolerance to accommodate for uncertainties in target data. With respect to these observables, different parameterizations with different "phenotypes" can have the same objective cost, as long as they are within the specified tolerances.

Bonded interactions: Table 4 lists the range of permitted bond parameters used in the optimization. The resulting bonded parameters of P g are shown in Figure 6. For equilibrium bond lengths b 0 there is little variation between different candidate solutions. This strong consensus suggests that the optimization has converged and that small changes in equilibrium bond length are linked to significant cost changes. The situation for the force constants is quite different. The values fluctuate over a relatively large range, compared to the predefined domain of permitted values. The measured bond length distributions (Figure 4A) show that these seemingly substantial differences in force constant values have only minor effects on the molecule's behavior.

The situation for the angle FF parameters is similar. The equilibrium values show smaller variances than the force constants, compared to their respective domain sizes of applicable values. Again, the differences in FF parameters have little effect on the observed distributions (cf. Figure 4B). Notably, the optimal force constants for the angles PO4-AM1-T1A and AM2-C1B-C2B were close to or at the maximum of their permitted ranges. Further optimization was therefore likely hindered, and a wider range should have been chosen.

In a similar vein to the discussion surrounding non-bonded parameters, the relatively wide range of force constants in P g indicates that additional metrics or training systems could be employed to further optimize the overall performance of candidate solutions while maintaining the quality of the employed observables. For instance, exploring lipids in environments other than a bilayer, which induce different lipid conformations, could benefit from a candidate solution with a lower angle force constant to allow for increased conformational variation. ). In the optimization runs with resampling an initial computational budget of 64 • 3 = 192 MD simulation slots was used for one full objective function evaluation of each particle, and a second equally sized computational budget was allocated to reevaluate the melting temperature (the target observable with the largest variance) of the best 16, best 32, or all 64 candidate solutions of the current iteration. For brevity we will refer to these as mv-PSO-R16, mv-PSO-R32, and mv-PSO-R64. Due to the fixed computational budget, for each particle involved in resampling, T m was reevaluated 12, 6, or 3 times. As half of the total computational budget was used for resampling, the number of iterations was set to 42 in these runs.

From the literature on PSO noise-mitigation [START_REF] Rada-Vilela | Population statistics for particle swarm optimization: Resampling methods in noisy optimization problems[END_REF][START_REF] Rada-Vilela | Deception, blindness and disorientation in particle swarm optimization applied to noisy problems[END_REF] we draw the expectation that which of the resampling, or no resampling, strategies is the best, depends on the level of noise. If noise levels are very low, the additional number of possible iterations, when forgoing resampling, could lead to better solutions. For intermediate noise levels, initial fitness evaluation results in a sufficient differentiation of good and bad solutions, i.e., overall sorting is roughly correct, and the focus on improving sorting of the very best solutions is most helpful. In case of even higher noise levels initial sorting would be vastly incorrect and a larger fraction of the swarm needs to be resampled to achieve satisfactory overall sorting. As a consequence, the sorting quality of the very top would be degraded, as there is less computational budget allocated here.

The true quality of a candidate solution is not necessarily reflected by the cost estimated during an optimization run, as there is some uncertainty in estimates of target observables other than T m , and the confidence level of the T m estimation with different resampling settings differs vastly. Therefore, validation is required. As we are mostly interested in the quality verification of the best solutions, the first step of the screen-to-the-best procedure from Boesel et al. [START_REF] Boesel | Using Ranking and Selection to "Clean up" after Simulation Optimization[END_REF] can be used to select the statistically equivalent set of candidate solutions. For mv-PSO-R16 the set P g t contains 69 candidate solutions. Due to the increased uncertainty in mv-PSO-R32 and mv-PSO-R64, their respective sets P g t contain hundreds of candidate solutions. To keep the computational cost for validation manageable, we selected only the 72 best solutions of these optimization runs for validation. As there are no variance estimates in the optimization run without resampling, the selection procedure is not applicable. Again, the 72 best solutions from the optimization run were selected for validation. All candidate solutions chosen for validation were fully (all training systems, all observables) reevaluated 20 times. The resulting rerun cost vs. the originally estimated cost is shown in Figure 7.

Clearly, mv-PSO-R16 gave the best results, while the quality of the best solutions in the three other cases does not differ much. Furthermore, the fact that for all selected candidate solutions of mv-PSO-R0 the rerun cost estimate is substantially higher than the original cost estimate indicates that these original estimate are particularly favorable. While there are also candidate solutions with substantial differences in original and rerun cost for the resampling systems -in this case mostly caused by APL fluctuations -these are much less frequent and there is much better correlation between original and rerun cost (Pearson correlation coefficient 0.21 vs. 0.64, for mv-PSO-R0 and mv-PSO-R16, respectively).

Our interpretation of these results is the following: The noise level is low enough, so that even without noise-mitigation, the sorting of candidate solutions is correct in a coarser sense and the swarm is guided towards the "correct" vicinity in parameter space. Yet, noise levels are substantial enough, so that resolution of finer cost differences is impeded. Only the concentrated allocation of the resampling budget on the top 16 solutions lowers the cost original cost We have illustrated how to apply mixed-variable particle swarm optimization for automated CG molecule parameterization. As an example application, we parameterized the sphingolipid linker region for the Martini 3 FF. The newly parameterized sphingomyelin model reproduces important target observables accurately, including the melting temperature, which was ≈ 20 K off target before and is now within ≈ 2 K of the experimental reference. Notably, reproduction of experimental melting temperatures had been historically problematic in Martini lipid models. [START_REF] Marrink | Simulation of gel phase formation and melting in lipid bilayers using a coarse grained model[END_REF] The mixed-variable approach offers a major advantage when parameterizing molecules for building-block force fields. Due the explicit use of building blocks, every candidate model is a valid parameterization in the given FF. Otherwise, changing non-bonded interaction parameters of the FF's building blocks breaks the validity of their parameterization. Candidate solutions generated by a continuous treatment of non-bonded interactions have to be converted to a valid FF model, followed by additional validation of this model.

A drawback of the mixed-variable treatment is that some advanced improvements to PSO, such as fuzzy parameter tuning of Nobile et al., [START_REF] Nobile | Fuzzy Self-Tuning PSO: A settings-free algorithm for global optimization[END_REF] are not directly applicable to mv-PSO, because in the categorical representation there is no similarity metric, which is utilized in the PSO parameter tuning. This could be overcome by using discrete ordered representation for non-bonded interactions instead of the categorical treatment.

One of the great benefits of automated parameterization algorithms is the simultaneous optimization against multiple structural and thermodynamic target data. As thermodynamic observables can be expensive to estimate accurately in MD simulations, the formal consideration of noise in objective function values is an important conceptual improvement.

As demonstrated, optimization with applied noise-mitigation produced significantly better solutions and the utilized screen-to-the-best procedure provides a systematic approach to the post-optimization selection of the best model.

Although we have demonstrated the adverse effects of objective function value noise on the sorting and performance of PSO, it is important to note that the non-deterministic nature of particle swarm optimization necessitates multiple repetitions of full optimization runs to confidently determine the most effective noise-mitigation setting. Achieving a high level of confidence in identifying the optimal approach would require a significant number of iterations. Furthermore, the 'ground truth', i.e., the true score of a candidate parameterization, is unknown, hence a large amount of validation simulations would be required. This is not feasible, due to a high computational cost. Rigorous development and testing of noise-mitigation strategies should not be done with objective function evaluations that require costly MD simulations, and are therefore beyond the scope of this paper. Moreover, the additionally gained insight, would only be of moderate value. The PSO literature has shown that under significant noise PSO performance is degraded and performance differences between resampling methods for noise-mitigation are problem and noise-level dependent. Generally, noise-mitigation methods employing OCBA perform the best under various circumstances, [START_REF] Rada-Vilela | Population statistics for particle swarm optimization: Resampling methods in noisy optimization problems[END_REF][START_REF] Rada-Vilela | Population statistics for particle swarm optimization: Hybrid methods in noisy optimization problems[END_REF] but its sequential secondary budget allocation puts constraints on the parallelization of the parameterization algorithm. Still, its integration into the parameterization pipeline should be explored in the future.

Together with the general benefits of automation, the here-presented conceptual advantages will further facilitate rigorous CG molecule parameterization. The CGCompiler Python package that comes with our method is tailor-made for parameterization tasks in buildingblock FFs, such as Martini. Also larger building blocks, i.e., a molecule class with shared regions can be parameterized simultaneously. Our approach is not limited to lipid parameterization, but can be applied to any kind of molecule. CGCompiler can be easily adapted to the needs of a specific parameterization task. Implementing new observables is not much different from writing Python functions for analyzing MD data. Importantly, our automation platform eases collaborations between individual researchers since a clear overview of the parameterization flow is provided. This also renders force-field reproducibility as well as retrospective force-field corrections, such as corrections to the targets (e.g., improved atomistic force-fields or simulation settings) or inclusion of additional targets rather straightforward.

The here-presented study focuses on method development and the sphingolipid linker parameterization was merely a test case. The parameters of the head group and lipid tails, predefined in our study, are still actively improved/(re)parameterized by the core developers. [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] Once these final parameters are released, reparameterization of the linker may be necessary, ideally with an even broader set of training systems, including liquid ordereddisordered phase behavior.

In order to achieve fully automated molecule parameterization in high-throughput applications, the development of an automated mapping and selection of bonded terms remains a crucial component. Currently, mapping and parameter optimization are separate tasks, but integrating an automated mapping scheme into the parameterization pipeline could be facilitated prior to employing mixed-variable particle swarm optimization, utilizing CGCompiler.

The choice of bonded parameters not only influences the accuracy of the model but also impacts simulation stability. Various strategies, such as the use of virtual sites, restricted bending potentials, hinge and "divide and conquer" constructions, [START_REF] Bulacu | Improved Angle Potentials for Coarse-Grained Molecular Dynamics Simulations[END_REF][START_REF] Alessandri | Martini 3 Coarse-Grained Force Field: Small Molecules[END_REF] have been previously described to address instability. Additionally, careful consideration of constraints is necessary to ensure simulation stability and prevent artificial temperature gradients. [START_REF] Fábián | Optimal Bond Constraint Topology for Molec[END_REF][START_REF] Thallmair | Nonconverged Constraints Cause Artificial Temperature Gradients in Lipid Bilayer Simulations[END_REF] These aspects should be incorporated as essential steps in a future fully automated parameterization pipeline.

Another future prospect is the advancement of true non-scalarized multi-objective optimization, which eliminates the need for user-defined weights on the targets within the objective function. However, it can also be argued that these user-defined weights, which reflect the importance of targets based on intuition, experience, or additional knowledge, along with the predefined set of relevant structural and thermodynamic targets for the CG force-field, encompass what is commonly known as the "force-field's philosophy". In this sense, the user-defined weights embody the guiding principles that shape the force-field.

  was used for electrostatic interactions. Temperature coupling was done via the velocity rescale algo-16 https://doi.org/10.26434/chemrxiv-2023-r97bm ORCID: https://orcid.org/0000-0001-5239-7124 Content not peer-reviewed by ChemRxiv. License: CC BY-NC 4.0 rithm 45 with a coupling time τ t = 1.0 ps. System pressures were held at 1 bar by using the Parinello-Rahman barostat[START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] with a coupling time τ p = 5.0 ps. Pressure coupling was applied isotropically for aqueous solutions and semi-isotropically for membrane systems.

Figure 3 :

 3 Figure 3: Thickness, average area per lipid and melting temperature for the set of statistically equal candidate solutions that remained after the second screen-to-the-best procedure performed after reevaluating the initial set 20 times.

Figure 4 :Figure 5 :

 45 Figure 4: Validation of targets from rerun simulations for the set P g . A) Bond length distributions. B) Angle distributions.

Figure 6 :

 6 Figure 6: Force field parameters of the set of statistically equivalent solutions P g for the sphingolipid linker region. A) Bead probability distributions. B) Bond parameters. Dashed lines are upper and lower parameter limits. C) Angle parameters. Dashed lines are upper and lower parameter limits. The equilibrium angles of AM1-T1A-C2A and AM2-C1B-C2B are not varied during optimization. They are fixed at 180 • .

Figure 7 :

 7 Figure 7: Comparison of cost estimated during the optimization run and average cost estimated from repeated reruns of P g t in (A) and the 72 best candidate solutions in (B)-(D). Error bars are standard errors. (A) Original cost: 16 particles resampled, 1+12 T m samples. (B) Original cost: 32 particles resampled, 1+6 T m samples. (C) Original cost: 64 particles resampled, 1+3 T m samples. (D) No resampling during optimization, but twice as many iterations

  1 and t β,ν is the β quantile of the t distribution with ν degrees of freedom 2. Set P g t = {i :i ∈ P t , fi ≤ fj + W ij , ∀ i = j ∈ P t } 3. Return P g t W ijis the half-width of pooled t-confidence intervals on the difference between the scores of solutions i and j.19 Therefore, the procedure entails a pair-wise comparison of solutions and determines if differences of the sample averaged scores are statistically significant.19

zation

Figure 2: CG description of sphingomyelin and ceramide.

Table 1 :

 1 Atomistic target system details. In the naming scheme of the CHARMM FF, SSM and CHL1 denote sphingomyelin (18:0) and cholesterol, respectively.

	system	lipids	# TIP3P # NA # CL T / K sim. time / ns
	DPSM128 328K	128 SSM	5120	-	-	328.15	150
	POPC SSM CHOL 100 POPC 100 SSM 100 CHL1	9000	18	18	321.15	300

Table 2 :

 2 Coarse-grained training system details.

	system	lipids	# W # NA # CL	T / K
	DPSM128 328K	128 DPSM	1177	-	-	328.15
	DPSM256 biphasic 256 DPSM half gel/half liquid 2300	26	26	286, 291, 296, 301, 303,
						305, 307, 308, 309, 310,
						311, 316, 321, 326
	POPC SSM CHOL 96 POPC 96 DPSM 96 CHOL 2124	23	23	321.15

Table 3 :

 3 Weights of observables w o and system specific observable weights w o,s for optimization run 1.Figure3shows thickness, average area per lipid and melting temperature of pure DPSM membranes for the set of statistically equal candidate solutions that remained after the second screen-to-the-best procedure performed after reevaluating the initial set 20 times. All new candidate solutions outperform the current DPSM model regarding thickness. The average area per lipid of the current model is closer to the target value, but most of the candidate solutions are within the tolerance of 1.5% deviation. In general, thickness and APL are inversely correlated, increasing one will always result in decreasing the other, there-

	observable	w o	w o, DPSM128 w o, DPSM256 w o, POPC SSM CHOL
	bond length dist.	1	1	0	1
	angle dist.	100	1	0	1
	d HH	500	1	0	0.25
	APL	1000	1	0	0.25
	T m	250	0	1	0
	RDF COM DPSM-CHOL	1	0	0	1

fore, with both values inside the tolerance, the new models represent a better balance of thickness and APL. It is important to note that in the comparison, SM(18:0) was used as the atomistic target. The current tail model of the Martini FF represents both SM(16:0) and SM(18:0). The CHARMM model for SM(16:0) exhibits a reduced thickness when compared

Table 4 :

 4 Bonded interactions. GROMACS function type; permitted parameter ranges for equilibrium bond length / angle, and corresponding force constants.

		GROMACS
	bond	bond func. type	b 0 / nm	fc / kJ/mol/nm 2
	PO4-AM1	1	0.25 -0.40	1000 -9000
	AM1-AM2	1	0.20 -0.35	1000 -9000
	AM1-T1A	1	0.40 -0.55	1000 -9000
	AM2-C1B	1	0.25 -0.50	1000 -9000
		GROMACS
	angle	angle func. type	a 0 / deg	fc / kJ/mol
	PO4-AM1-AM2	2	90 -180	5 -100
	PO4-AM1-T1A	2	90 -180	5 -100
	AM1-T1A-C2A	2	180	5 -100
	AM2-C1B-C2B	2	180	5 -100
		25

https://doi.org/10.26434/chemrxiv-2023-r97bm ORCID: https://orcid.org/0000-0001-5239-7124 Content not peer-reviewed by ChemRxiv. License: CC BY-NC 4.0

https://doi.org/10.26434/chemrxiv-2023-r97bm ORCID: https://orcid.org/0000-0001-5239-7124 Content not peer-reviewed by ChemRxiv. License: CC BY-NC 4.0

https://doi.org/10.26434/chemrxiv-2023-r97bm ORCID: https://orcid.org/0000-0001-5239-7124 Content not peer-reviewed by ChemRxiv. License: CC BY-NC 4.0 estimation errors sufficiently such that improved candidate solutions can be found.

Acknowledgement K.S.S. and H.J.R. thank the NWO Vidi Scheme, The Netherlands, (project number: 723.016.005) for funding this work. H.J.R thanks the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for funding this work under Germany's Excellence Strategy -EXC 2033 -390677874 -RESOLV. K.S.S. and H.J.R. gratefully acknowledge the computing time granted by the Resource Allocation Board and provided on the supercomputer Lise and Emmy at NHR@ZIB and NHR@Göttingen as part of the NHR infrastructure. The calculations for this research were conducted with computing resources under the project nip00054. K.S.S. and H.J.R. gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). P.C.T.S. acknowledges the support of the French National Center for Scientific Research (CNRS) and the funding from a research collaboration agreement with PharmCADD. L.M. acknowledges funding by the Institut National de la Santé et de la Recherche Médicale (INSERM)

Data availability

The CGCompiler Python package is available upon request from the authors and will be made publicly available on github at a later stage.

A Gromacs topology file of the final sphingomyelin (and ceramide) can be download from the Martini Database server 63 (https://mad.ibcp.fr).

Supporting Information Available

The following files are available free of charge.

• SI.pdf: Supporting Information including methods to estimate melting temperature additional plots regarding noise -Sphingomyelin-cholesterol 2d-center-of-mass radial distribution functions -DPSM Gromacs topology file


Eugenie Huillet
email: eugenie.huillet@jouy.inra.fr

Michel Gohar
email: michel.gohar@jouy.inra.fr

Chapter 21.5. Quorum sensing in Bacillus cereus in relation with cysteine metabolism and oxidative stress response

Keywords: Bacillus cereus, quorum sensing, AIP, PlcRa, cysteine metabolism, CymR, AI-2, Lsr-like system, oxidative stress response, biofilm, global transcriptional analysis

A number of peptide-mediated quorum sensing (QS) systems regulate important biological functions in bacteria of the Bacillus cereus group. Sporulation, virulence and biofilm formation are controlled by such systems in these bacteria. PapR, a secreted heptapeptide, activates PlcR, the main transcriptional regulator of virulence in B. thuringiensis and B. cereus. We recently found that PlcRa, a PlcR paralogue, also works in QS and is involved in oxidative stress response and cysteine metabolism during the early stationary phase of bacterial growth. Notably, we have reported that PlcRa activates the expression of several genes involved in the pathway of cysteine synthesis from methionine including the S-adenosylmethionine (SAM) recycling pathway. This transcriptional control of SAM recycling pathway by PlcRa impacts another QS, the autoinducer-2 (AI-2) system. AI-2 is a signal metabolic molecule produced by LuxS, a central metabolic enzyme involved in SAM recycling, found in many bacterial species and thus proposed to enable interspecies communication. We have reported that B. cereus synthetizes and recognizes AI-2 as an extracellular signal. Most particularly, we have shown that AI-2 inhibits biofilm formation in a concentration-dependent manner.

In this review, we provide an overview of these two QS systems in B.cereus species (signal production, detection, transduction) with the description of AI-2 connection. We highlight how CymR regulon -involved in the global regulation of gene expression in response to cysteine availability-together with PlcRa system is important for AI-2 production and present the role of PlcRa in oxidative stress response.

Introduction

Quorum sensing (QS) allows bacteria to share information about cell density to adjust gene expression accordingly. This cellcell communication relies on secreted signal molecules, called autoinducers (AI) or pheromones. This system enables bacteria to undergo energetically expensive processes as a collective only when the impact of these processes on their environment is more likely to be efficient [START_REF] Dunny | Cell-cell communication in grampositive bacteria[END_REF].

Intercellular communication processes are crucial in coordinating growth, adaptation, pathogenicity, biofilm formation and sporulation in B. cereus group [START_REF] Slamti | Quorum sensing in Bacillus thuringiensis is required for completion of a full infectious cycle in the insect[END_REF]. The B. cereus group includes well-known spore-forming pathogens of mammals (B. anthracis and B. cereus) and insects (B. thuringiensis). B. cereus is frequently associated with foodborne infections causing gastroenteritis [START_REF] Stenfors Arnesen | From soil to gut: Bacillus cereus and its food poisoning toxins[END_REF]. B. cereus has also been described as the cause of nosocomial infections in immuno-compromised patients and is involved in endocarditis, endophtalmites or meningitis that have proved fatal [START_REF] Stenfors Arnesen | From soil to gut: Bacillus cereus and its food poisoning toxins[END_REF]. The capacity of B. cereus to sporulate allows this bacterium to resist the usual cleaning procedures used in the food industry, resulting in the presence of B. cereus in many raw and processed foods, such as rice, spices, milk, vegetables, meats and various desserts [START_REF] Abee | Germination and outgrowth of spores of Bacillus cereus group members: diversity and role of germinant receptors[END_REF].

The PlcR/PapR QS system controls the expression of genes encoding exported virulence factors, including degradative enzymes, enterotoxins and hemolysins in B.cereus and B.thuringiensis species [START_REF] Agaisse | PlcR is a pleiotropic regulator of extracellular virulence factor gene expression in Bacillus thuringiensis[END_REF]. The transcription of these virulence genes in bacteria growing in a rich medium is activated at the onset of stationary phase or transition phase by the active quorum sensor PlcR. PlcR is activated by binding to PapR, a signaling peptide produced as a propeptide under the control of PlcR (Figure 21.5.1). PapR undergoes extracellular processing, to generate an active heptapeptide [START_REF] Bouillaut | Molecular basis for group-specific activation of the virulence regulator PlcR by PapR heptapeptides[END_REF], which is then re-imported into the bacterial cell via the oligopeptide permease system, OppABCDF. Within the cell, PapR interacts with PlcR and the resulting complex binds PlcR target sites on DNA [START_REF] Slamti | A cell-cell signaling peptide activates the PlcR virulence regulon in bacteria of the Bacillus cereus group[END_REF], resulting in the activation of the PlcR regulon, which contains 45 genes [START_REF] Gohar | The PlcR virulence regulon of Bacillus cereus[END_REF]. For a more detailed overview of PlcR/PapRa system we recommend a recent review by Slamti [START_REF] Slamti | Quorum sensing in Bacillus thuringiensis is required for completion of a full infectious cycle in the insect[END_REF] Determination of the structure of PlcR resulted in the identification of a new family of central regulatory quorum sensors (the RNPP family) found exclusively in Gram-positive bacteria with a low G+C content [START_REF] Declerck | Structure of PlcR: Insights into virulence regulation and evolution of quorum sensing in Gram-positive bacteria[END_REF]. These quorum sensors include PlcRa and NprR from B. cereus [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF][START_REF] Perchat | A cell-cell communication system regulates protease production during sporulation in bacteria of the Bacillus cereus group[END_REF], PrgX from Enteroccocus faecalis [START_REF] Shi | Structure of peptide sex pheromone receptor PrgX and PrgX/pheromone complexes and regulation of conjugation in Enterococcus faecalis[END_REF] and RAP phosphatases from Bacillus subtilis [START_REF] Parashar | Structural basis of response regulator dephosphorylation by Rap phosphatases[END_REF][START_REF] Solomon | Purification and characterization of an extracellular peptide factor that affects two different developmental pathways in Bacillus subtilis[END_REF]. All these RNPP regulators are activated through a secreted signaling peptide that interacts with the TPR activation domain.

We recently characterized PlcRa, which activates the oxidative stress response and cysteine metabolism in transition state cells [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. PlcRa displays about 29% sequence identity and about 50% similarity to PlcR. We found that the PlcRa protein is a 3D structural homolog of PlcR, suggesting that the TPR peptide-binding domain may be arranged similarly to form a pocket that is responsible for peptide binding. PlcRa is activated by the product of the papRa gene, a secreted signaling peptide (Figure 21.5.1) [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. Biochemical and genetic analyses demonstrated the role of the heptapeptide PapRa 7 in PlcRa activity. Purified PlcRa specifically binds to a PlcRacontrolled promoter in the presence of PapRa heptapeptide (abrB2 gene) [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. The production, the maturation of However, a conserved quorum-sensing system is found both in Grampositive and Gram-negative bacteria [START_REF] Pereira | AI-2-mediated signalling in bacteria[END_REF]. The signal molecule of this system is AI-2 (autoinducer 2), a detoxification product of the activated methyl cycle or SAM recycling linked to cysteine biosynthesis from methionine (Figure 21.5.2). AI-2 is a derivative of (S) 4,5-dihydro-2,3-pentanedione (DPD). DPD itself is a by-product of the activated methyl cycle (AMC) which generates the methyl donor, S-adenosylmethionine (SAM), required for a number of metabolic reactions including Once in the cell, AI-2 is phosphorylated by the LrsK kinase, enabling its binding to, and the inactivation of, the LrsR transcriptional repressor. Genes regulated by LrsB can therefore be transcribed (Figures 21.5.3,21.5.4).

Despite the presence of PfS and LuxS in most bacteria, not all of these bacteria use AI-2 as a signal molecule for quorum-sensing purposes, since some of them lack the AI-2 receptor [START_REF] Rezzonico | Lack of genomic evidence of AI-2 receptors suggests a non-quorum sensing role for luxS in most bacteria[END_REF]. However, pathogens of the B. cereus group (B. anthracis, B. cereus and B. thuringiensis) possess the genetic equipment for AI-2 secretion, importation and processing. In addition, the LuxS-dependent production of AI-2 and its effect on growth has been demonstrated in B. anthracis [START_REF] Jones | Detection of a luxS-signaling molecule in Bacillus anthracis[END_REF]. Cysteine metabolism is well characterized in B. subtilis [START_REF] Burguiere | Three different systems participate in L-cystine uptake in Bacillus subtilis[END_REF][START_REF] Even | Global control of cysteine metabolism by CymR in Bacillus subtilis[END_REF][START_REF] Hullo | Conversion of methionine to cysteine in Bacillus subtilis and its regulation[END_REF].

In this bacterium, cysteine is synthesized either from homocysteine via the reverse transsulphuration pathway [START_REF] Hullo | Conversion of methionine to cysteine in Bacillus subtilis and its regulation[END_REF] In B. subtilis [START_REF] Even | Global control of cysteine metabolism by CymR in Bacillus subtilis[END_REF][START_REF] Tanous | The CymR regulator in complex with the enzyme CysK controls cysteine metabolism in Bacillus subtilis[END_REF] CymR for "cysteine metabolism repressor" is the master regulator of cysteine metabolism. CymR represses the transcription of a large set of genes involved in cystine uptake and cysteine biosynthesis. The CymR repressor belongs to the widespread and poorly characterized Rrf2 family of transcription factors [START_REF] Tanous | The CymR regulator in complex with the enzyme CysK controls cysteine metabolism in Bacillus subtilis[END_REF][START_REF] Shepard | Insights into the Rrf2 repressor family--the structure of CymR, the global cysteine regulator of Bacillus subtilis[END_REF]. CymR displays a particular mechanism of control of its activity, with the formation of a regulatory complex with CysK, an OAS-thiol-lyase. In this complex, CymR is the DNA-binding protein, and CysK increases the stability of the interaction between CymR and DNA [START_REF] Tanous | The CymR regulator in complex with the enzyme CysK controls cysteine metabolism in Bacillus subtilis[END_REF].

In the signal transduction pathway controlling cysteine metabolism, CysK, via its substrate OAS, is the sensor of the cysteine pool in the cell for the regulatory complex, and transmits this information to the CymR repressor [START_REF] Hullo | Conversion of methionine to cysteine in Bacillus subtilis and its regulation[END_REF]. OAS is produced from acetyl-CoA and serine by CysE, the serine acetyltransferase while the OAS-thiol-lyase, CysK, further condenses sulphide and OAS to form cysteine. Cysteine limits the OAS pool via its negative control of CysE synthesis and activity.

In the presence of cysteine, the intracellular concentration of OAS is low, the CymR-CysK complex is formed, and transcriptional repression of the CymR regulon occurs (Figure 21.5.6.A). Under conditions of cysteine limitation, the pool of OAS increases until it reaches a concentration preventing the formation of the CymR-CysK complex or dissociating the complex in vitro (Figure 21.5.6.B). Thus, no repression of transcription of the CymR regulon occurs when cysteine is limiting [START_REF] Even | Global control of cysteine metabolism by CymR in Bacillus subtilis[END_REF][START_REF] Hullo | Conversion of methionine to cysteine in Bacillus subtilis and its regulation[END_REF][START_REF] Tanous | The CymR regulator in complex with the enzyme CysK controls cysteine metabolism in Bacillus subtilis[END_REF]. PlcRa controls the expression of CymR regulon. We characterized the regulatory role of PlcRa, by comparing the transcriptomes of the ∆plcRa and wild-type B. cereus ATCC 14579 strains during early stationary phase in a complex rich medium, the Luria Bertani broth. Assessments were carried out one hour after the onset of stationary phase (t 1 ), when plcRa expression increases, and two hours after the onset of stationary phase (t 2 ), when plcRa expression has reached a plateau [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF].

Twenty of the 49 genes upregulated by PlcRa encode proteins involved in sulfur metabolism (Table 21.5.1). These genes were found to be differentially expressed two hours after the onset of stationary phase [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. The sulfur metabolism genes of B. cereus remain poorly annotated, despite the availability of several B. cereus genomes. We therefore reconstructed the sulfur metabolism pathway, by searching for orthologs of B. subtilis genes in the B. cereus ATCC 14579 genome [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF] (Figure 21.5.5). The transport of sulfur sources [START_REF] Burguiere | Three different systems participate in L-cystine uptake in Bacillus subtilis[END_REF] and the two major cysteine biosynthetic pathways in B. subtilis -the thiolation pathway, which requires sulfide, and the reverse transsulfuration pathway, which converts homocysteine to cysteine, with cystathionine formed as an intermediate [START_REF] Even | Global control of cysteine metabolism by CymR in Bacillus subtilis[END_REF][START_REF] Hullo | Conversion of methionine to cysteine in Bacillus subtilis and its regulation[END_REF] -are conserved in B. cereus ATCC 14579. We were able to identify all the enzymes and transporters required for these pathways other than those for the reduction of sulfite to sulfide [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF] No repression of transcription of the CymR regulon occurs when cysteine is limiting in B.subtilis. We hypothesized that PlcRa activates CymR-controlled genes when cysteine is limiting during the transition phase in B.cereus. Indeed, we found that the expression of the PyrrT'-lacZ fusion was induced at the onset of stationary phase in the wild-type strain, suggesting cysteine depletion in the growth medium leading to the induction of the methionine-to-cysteine conversion pathway [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF].

No such induction was observed in the plcRa mutant strain [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]

Role of the cysteine metabolism in oxidant protection in

Bacillus. Many agents that restrict bacterial growth and developmental behavior induce the production and/or accumulation of reactive substances, such as reactive oxygen species (ROS), reactive nitrogen species (RNS), and toxic electrophiles.

Increased levels of reactive species accompany stress, in many of its manifestations, and lead to macromolecular damage, an increased rate of mutagenesis, and cell death (Zuber,2009).

Bacillus respond to oxidative stress by the activation of different cellular defence mechanisms [START_REF] Zuber | Management of oxidative stress in Bacillus[END_REF]. These are composed of scavenging enzymes, as well as protection and repair systems presumably organized in highly sophisticated networks [START_REF] Zuber | Management of oxidative stress in Bacillus[END_REF][START_REF] Mols | Primary and secondary oxidative stress in Bacillus[END_REF]. Cysteine alone, cysteine molecules and cysteine proteins are major actors in oxidant protection. In particular, cysteine is the direct precursor of the low-molecular mass antioxidant molecule, Bacillithiol, a key actor of the maintenance of the cytosolic redox balance and of the adaptation to reactive oxygen species in Bacilli [START_REF] Newton | Bacillithiol is an antioxidant thiol produced in Bacilli[END_REF][START_REF] Helmann | Bacillithiol, a new player in bacterial redox homeostasis[END_REF]. In addition, S-cysteinylation is a general mechanism for thiol protection of B. subtilis proteins after oxidative stress [START_REF] Hochgrafe | S-cysteinylation is a general mechanism for thiol protection of Bacillus subtilis proteins after oxidative stress[END_REF]. Consequently, the CymR regulon involved in cysteine biosynthesis is involved in oxidant protection. It has been shown that a cymR mutant is highly sensitive to H 2 O 2 , disulfide, paraquat, copper and tellurite induced stresses in B. subtilis [START_REF] Soutourina | CymR, the master regulator of cysteine metabolism in Staphylococcus aureus, controls host sulphur source utilization and plays a role in biofilm formation[END_REF]. The CymR regulon is induced during disulfide or superoxide stresses in B. subtilis [START_REF] Mostertz | Transcriptome and proteome analysis of Bacillus subtilis gene expression in response to superoxide and peroxide stress[END_REF].

PlcRa: role in the oxidative stress response.

PlcRa is mainly involved in the activation of genes that function in oxidant protection through the control of the expression of CymR regulon [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. PlcRa also activates the expression of ahpCF operon encoding the alkyl hydroperoxyreductase (AhpR), a detoxification system composed of two enzymes that catalyses the breakdown of molecular oxygen to hydrogen peroxide (H 2 O 2 ) which is reduced by the second enzyme to water. PlcRa activates the expression of dps2 gene encoding a Dps-like mini-ferritin protein (BC5044) which is homologous to the Dps2 protein of B.anthracis that has recently been shown to play a major role in oxidative stress resistance [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. These genes are induced during H 2 O 2 stress in B. cereus [START_REF] Ceragioli | Comparative transcriptomic and phenotypic analysis of the responses of Bacillus cereus to various disinfectant treatments[END_REF].

We tested the effect of plcRa inactivation on the susceptibility of B. cereus to H 2 O 2 and diamide stresses. Diamide causes thiol oxidation and disulfide stress. The viability of the ∆plcRa mutant and the parental strain was tested after the addition of H 2 O 2 (1 mM). A 7-fold reduction in survival was observed for the ∆plcRa mutant compared to the wild-type strain. The complementation of the ∆plcRa mutant restored a wild-type phenotype. Then, the viability of the ∆plcRa mutant and the parental strain was tested after the addition of diamide (10 mM).

A 4-fold reduction in survival was observed for the plcRa mutant compared to the wild-type strain while the introduction of plcRa in trans restored a wild-type phenotype [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF].

These results showed that plcRa inactivation led to an increased sensitivity to H 2 O 2 and disulfide stresses. We also investigated the effects of cystine addition during peroxide stress [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. Cystine (1 mM) was added to the culture in midexponential growth phase (figure 21.5.8). Two hours later, at the end of exponential phase, H 2 O 2 (0.4 mM) was added to the medium.

Both strains presented a growth arrest one hour after the H 2 O 2 addition which was characterized by a OD 600 measurements drop. In addition, in the presence of H 2 O 2 , OD 600 measures of the ∆plcRa mutant strain were lower than the wild type strain. When cystine was added, the growth arrest for both strains was abolished. Thus, cystine significantly reduced the sensitivity of these cells to H 2 O 2 stress in our growth conditions. Moreover, these results strongly suggest that cystine transport is efficient in the ΔplcRa mutant, as in the wild type and it might be due at least partly to the TcyABC system (BC0872-BC0873-BC0874) that is not controlled by PlcRa (figures 21.5.5, 21.5.6). There is accumulating evidence suggesting that metabolic imbalances also induce bacterial processes that are normally associated with the oxidative stress response and/or with the general stress response. Altered transcriptional activites in plcRa mutant (117 genes differentially expressed with a more than three-fold difference) probably resulted in metabolic changes which seem to have caused general and oxidative stresses [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. Indeed, in the plcRa mutant, in LB medium without treatment, 10 genes that belong to the σ B regulon, are strongly higher expressed [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF]. This includes σ B and its associated regulatory proteins RbsV, RbsW and RbsP. In B.subtilis, elevated activity of σ B factor and induction of the σ B general stress regulon confer increased oxidant resistance [START_REF] Zuber | Management of oxidative stress in Bacillus[END_REF]. Upon exposure of B.cereus to mild and lethal H 2 O 2 concentrations the expression of numerous genes was affected, including genes involved in the σ B common response to general stresses [START_REF] Ceragioli | Comparative transcriptomic and phenotypic analysis of the responses of Bacillus cereus to various disinfectant treatments[END_REF]. Oxidative and disulfide stresses regulation of prophage genes has been reported in Lactobacillus plantarum, L. reuteri and B.subtilis [START_REF] Wilson | Transcriptional and metabolomic consequences of LuxS inactivation reveal a metabolic rather than quorum-sensing role for LuxS in Lactobacillus reuteri 100-23[END_REF]. Notably in a plcRa mutant, the expression of 31 prophage genes is higher [START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF].

AI-2 mediated signalling in B.cereus: Lsr-like system

The production of AI-2 has been shown in B. anthracis to be LuxSdependent, and this production is maximal a the end of the exponential phase of growth and decreased in stationary phase [START_REF] Jones | Detection of a luxS-signaling molecule in Bacillus anthracis[END_REF]. In B. cereus too, production of AI-2 increased in exponential phase and decreased in stationary phase (Figure 21.5.8), suggesting that this quorum-sensing system is active during the transition between the exponential and the stationary phases of growth. In the same species, the effect of exogenously added AI-2 in a microtiter plate biofilm assay was assayed. AI-2 inhibited in a dose-dependent way the production of biofilm in these conditions. It also induced a decrease in biofilm biomass when this compound was added to mature biofilms (Figure 21.5.9) [START_REF] Auger | Autoinducer 2 affects biofilm formation by Bacillus cereus[END_REF]. In Gram-positive as well as in Gram-negative bacteria, the AI-2 quorum-sensing system controls a number of phenotypes, including chemotaxis or toxin production in pathogens. However, the phenotype the most frequently targeted by this quorum-sensing system is biofilm formation [START_REF] Pereira | AI-2-mediated signalling in bacteria[END_REF].

Interestingly, while AI-2 promotes biofilm formation in most species, the inverse effect is observed in B. cereus. This effect might be of importance in the formation of multispecies biofilm, where AI-2, working as an interspecies communication molecule, could modulate the balance of biofilm species. While the lsr system appears to be fully functional in pathogenic bacteria of the B. cereus group, this is not a general rule in all bacteria of the Firmicutes phylum. In Listeria monocytogenes in Clostridium difficile or in Lactobacillus reuteri, for instance, the effect of luxS deletion is not a consequence of AI-2 absence, but is due to a perturbation of the bacterium metabolism (Challan [START_REF] Belval | Assessment of the roles of LuxS, S-ribosyl homocysteine, and autoinducer 2 in cell attachment during biofilm formation by Listeria monocytogenes EGD-e[END_REF][START_REF] Garmyn | Communication and autoinduction in the species Listeria monocytogenes: A central role for the agr system[END_REF][START_REF] Carter | Quorum sensing in Clostridium difficile: analysis of a luxS-type signalling system[END_REF].

Concluding remarks: interplay between PlcRa and Lsr-like

QS during the transition between vegetative and sporulation phases. At the end of the vegetative growth phase, bacterial cells face a number of challenges, including a decrease in the nutrient content of their environment and accumulation of oxidative products. Under these conditions, spore-forming bacteria may initiate sporulation, producing spores that can survive in unfavorable environmental conditions [START_REF] Piggot | Sporulation of Bacillus subtilis[END_REF]. Bacteria make use of various strategies to cope with environmental changes during the transition between the vegetative and sporulation phases [START_REF] Phillips | Bacillus subtilis sporulation and stationary phase gene expression[END_REF].

Alternatively to sporulation, bacteria can enter in another developmental state, the biofilm. The choice between sporulation and biofilm is likely to be a consequence of the speed at which Spo0A, the sporulation master regulator, is activated by kinases. Low Spo0A~P concentrations, while unable to initiate the sporulation cascade of events, can promote biofilm formation by repressing or activating the transcription of other regulator genes (namely abrB and sinI) [START_REF] Fujita | Evidence that entry into sporulation in Bacillus subtilis is governed by a gradual increase in the level and activity of the master regulator Spo0A[END_REF].

In this review, we presented a new peptide-mediated QS system, PlcRa/PapRa involved in AI-2 production and in the control of cysteine biosynthesis during the transition between vegetative and sporulation phases. The two QS systems, PlcRa and Lsr-like, are interconnected through AI-2 production. In this review, we highlight that cysteine metabolism is an essential compoment of Experiments were run in triplicates.

  PapRa and the binding of PapRa 7 with the TPR activation domain remain to established (Figure21.5.1).

Figure 21

 21 Figure 21.5.1 here...............................................

  methylations and polyamine biosynthesis. Release of methyl from SAM results in the production of S-adenosyl homocysteine (SAH), a toxic compound which conversion to DPD and homocysteine is catalysed by the PfS and LuxS enzymes. DPD spontaneously convert to AI-2 while homocysteine is recycled to SAM by the MetE, MetH and MetK enzymes. Therefore, AI-2 is both a product of the central metabolism and a signal molecule.

Figure 21

 21 Figure 21.5.2 here...............................................

Figure 21

 21 Figure 21.5.3 here...............................................

Figure 21

 21 Figure 21.5.4 here...............................................

  or from sulphide or thiosulphate via the thiolation pathway that directly incorporates these compounds into O-acetyl-L-serine (OAS) (Figure 21.5.5). Sulphide is obtained from the transport and reduction of inorganic sulphate or from organic sulphonates via their oxidation. Finally, methionine degradation is a biosynthesis pathway of cysteine formation and it proceeds via the formation of homocysteine using the S-adenosyl-methionine (SAM)recycling pathway followed by reverse transsulphuration, involving the MccA cystathionine β-synthase, and the MccB cystathionine γ-lyase (Figure 21.5.5).

Figure

  Figure 21.5.5 here………………………………………………………………………………………………….

Figure 21

 21 Figure 21.5.6.A here………………………………………………………………………………………………….

Figure 21

 21 Figure 21.5.6.B here………………………………………………………………………………………………….

A

  total of 18 PlcRa-controlled genes are involved in cysteine metabolism: transporters (tcyP) and biosynthesis of cysteine from sulfate (cysH operon and cysK gene) or methionine (yrrT-mtnN-mccAB operon and luxS gene) (Table 21.5.1 & Figure 21.5.5). BC4393, which shares 76 % identity with CymR from B. subtilis is probably the global regulator of cysteine metabolism in B.cereus. The CymR DNA-binding motif (WTAW2N4ATWNRWNANAN2GGAATT) is present in the promoter region of 7 PlcRa-controlled genes or operons involved in cysteine metabolism[START_REF] Huillet | PlcRa, a new quorum-sensing regulator from Bacillus cereus, plays a role in oxidative stress responses and cysteine metabolism in stationary phase[END_REF] (Figure21.5.7). A total of 14 PlcRa-controlled genes belongs to the CymR regulon (Table21.5.1).

Figure 21

 21 Figure 21.5.7 here…………………………………………………………………..

  . Here a proposed model for the global regulation of gene expression in response to cysteine availability in B.cereus (Figure 21.5.6AB). During the vegetative growth, in the presence of the preferred sulfur sources, sulfate or cysteine, CymR is active and represses the expression of CymRcontrolled genes. Upon sulfate or cysteine limitation, as presumably observed at the onset of stationary phase, CymR repressor is inactive, repression is thus relieved (see previous paragraph or legend Figure 21.5.6AB for details). In this growth condition, PlcRa activator, active during this transition phase, activates CymR-controlled genes.

Figure 21

 21 Figure 21.5.8 here...............................................

Figure 21

 21 Figure 21.5.9 here...............................................

B

  Figure Legends

Figure 21 .

 21 Figure 21.5.2. Biosynthesis of DPD. The transfert of a methyl

Figure 21 .

 21 Figure 21.5.3. Organisation of the lsr locus in B. cereus. The

Figure

  Figure 21.5.4. Schematic representation of two QS systems in B.

Figure 21 .

 21 Figure 21.5.5. Reconstruction of the sulfur metabolism pathway in

Figure

  Figure 21.5.6A-B. Schematic representation of the interplay

Figure

  Figure 21.5.7. Identification of a motif common to the promoter

Figure 21 .

 21 Figure 21.5.8. Cystine addition to a ∆plcRa culture improved

Figure 21

 21 Figure 21.5.9. Production and role of AI-2 in B. cereus (A)AI-2

  (Table 21.5.1 & Figure 21.5.5).

	Table 21.5.1 here……………………………………………..