Mariam Nakabuye

Abram Bunya Kamiza

Opeyemi Soremekun

Tafadzwa Machipisa

Emmanuel Cohen

Fraser Pirie

Oyekanmi Nashiru

Manjinder Elizabeth Young

S Sandhu

Ayesha A Motala

Manjinder S Sandhu

Tinashe Chikowore

Dr Segun Fatumo
email: segun.fatumo@lshtm.ac.uk

Genetic loci implicated in meta-analysis of body shape in Africans

Keywords: Body shape, Principal components, Anthropometric traits, Uganda, Zulu

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The prevalence of obesity has tripled since 1975, accounting for over 650 million obese individuals, globally. It is now one of the leading causes of non-communicable diseases, such as cardiovascular diseases (CVDs), metabolic diseases and cancer [START_REF] Oh | Effect of body shape on the development of cardiovascular disease in individuals with metabolically healthy obesity[END_REF][START_REF] Wang | Health and economic burden of the projected obesity trends in the USA and the UK[END_REF]. Anthropometric traits provide information about an individual's obesity through body shape, size and adiposity [START_REF] Amirabdollahian | Anthropometric Indicators of Adiposity Related to Body Weight and Body Shape as Cardiometabolic Risk Predictors in British Young Adults: Superiority of Waist-to-Height Ratio[END_REF]. Body shape indicates body fat distribution and different types of adiposity that can reveal an individual risk of developing certain diseases [START_REF] Jean-Pierre | Body Fat Distribution and Risk of Cardiovascular Disease[END_REF]. The gold standard methods of measuring adiposity include dual-energy X-ray absorptiometry (DXA) and magnetic resonance imaging (MRI) [START_REF] Borga | Advanced body composition assessment: from body mass index to body composition profiling[END_REF]. However, these methods are often not available in most low-income clinical settings and are also expensive to run. Whereas inexpensive methods of measuring body fat include the use of anthropometric traits such as weight, height, body mass index (BMI), waist circumference (WC), hip circumference (HC), waist to hip ratio (WHR) and a body shape index (ABSI), derived from weight, height, and WC.

Although inexpensive measures of body fat are readily available in most clinical settings from lower-to-high income and even rural-to-urban, these measures are often not robust enough to differentiate between fat and muscle mass, nor visceral adiposity and peripheral adiposity; therefore, poorly predicting body shape as an index for obesity. Thus, there is a growing interest to evaluate body shape derived from the principal component analysis (PCA) of anthropometric traits in large diverse datasets of uncorrelated variables that are combined to denote body shape. Only a few studies have applied a PCA approach to investigate genetic loci associated with a complex trait like body shape, predominately in Europeans [START_REF] Ried | A principal component meta-analysis on multiple anthropometric traits identifies novel loci for body shape[END_REF][START_REF] Aschard | Maximizing the Power of Principal-Component Analysis of Correlated Phenotypes in Genome-wide Association Studies[END_REF][START_REF] Fatumo | Complimentary Methods for Multivariate Genome-Wide Association Study Identify New Susceptibility Genes for Blood Cell Traits[END_REF]. However, due to the history of steatopygia and its elevated prevalence in Africans, a PCA approach assessing body shape as an index of obesity is needed in continental Africans.

To capture body shape in its multidimensional structure and to identify genetic loci associated with body shape in individuals of African ancestry, we applied PCA to several anthropometric traits, including weight, height, BMI, WC, HC and WHR. Herein, we performed J o u r n a l P r e -p r o o f a GWAS meta-analysis from two African cohorts to identify genetic loci associated with body shape as an index of obesity in individuals of African ancestry.

J o u r n a l P r e -p r o o f

Results

Overall, our study had 6,407 Ugandans from the general population cohort (GPC) and 2,598 South African Zulu cohort (SZC). The majority of our study participants were females in both the GPC (57.1%) and SZC (73.8%) datasets (Table 1). The median age of our study participants was 29.8 years in the Ugandan GPC and 52.0 years in the SZC, indicating SZC was predominately made up of older women. Moreover, we observed SZC had a higher median BMI, height, weight, WHR, HC, and WC compared with the Ugandans (Table 1).

PCA and body shape

In the Ugandans, our first PC explained 61.93% of the variation for all anthropometric traits. We noted that the first PC has high loadings for BMI, weight, HC and WC, thus capturing body fatness among Ugandans from GPC and we named it "weight driven body shape" (Fig. 1).

In the SZC, our first PC explained 56.41% of the variation and had high loadings for BMI, weight, HC and WC (Fig. 2). This PC also driven by weight, BMI , HC and WC, among the SZC and it was also named "weight driven body shape" to match with the one retained for the Ugandans. Based on the Kaiser rule, which chooses only PCs with eigenvalues >1 (to retain only PCs that explain more than one variable) [START_REF] Velicer | Component Analysis versus Common Factor Analysis: Some Further Observations[END_REF], we chose to retain the first and second PCs after PCA. We used the first PC of both Ugandans and Zulus, because of the similar anthropometric variables that were driving these principal components, which are the proxy for body shape in our analysis (Table S1, Fig S1). We then averaged loadings of first PC for Ugandans and Zulus to generate the PC which is representative of the most common body shape in Uganda and South African Zulu (Fig. S1). This first PC explained 60% of the variation for all anthropometric traits and had high loadings for weight, BMI, HC and WC; highlighting weight, BMI, HC and WC, as the main drivers of body shape in this context. The same approach was repeated for the sex stratified analysis. The male PC1 bodyshape was similar that of the combined analyisis, while the female pc1 body shape, had increseased loadings in the for waist-hip ration compared to the male and the combined PC1 body shape (Fig. 1) .

J o u r n a l P r e -p r o o f

Body shape association with cardiometabolic traits

We performed a linear regression to assess whether any changes in body shape in both the Ugandans and Zulus were associated with blood pressure (systolic and diastolic) and blood lipids (low-density lipoprotein cholesterol, triglycerides, total cholesterol, and highdensity lipoprotein cholesterol) traits in individuals of African ancestry. In our analysis, we observed that body shape changes were positively associated with blood pressure and blood lipids traits in Ugandans after adjusting for age and sex (Table.S2). Similarly, in the Zulus, the body shape was also positively associated with blood pressure traits, however, for blood lipid traits, the body shape was only positively associated with low-density lipoprotein cholesterol

(Table.S2).

Discovery of body shape genetic loci

We then performed genome-wide association analyses for body shape in Ugandans and Zulus and meta-analysed the summary statistics to identify genetic loci associated with body shape, as an index of obesity in both cohorts (Methods). Furthermore, we performed a meta-analysis using the summary statistics from the Ugandans and Zulus to identify additional genetic loci associated with body shape. We identified rs150717769 to be associated with body shape at genome-wide significance, P-value = 4.01e-09 (Table 2). The Manhattan and QQ plot for the meta-analysis of body shape in individuals of African ancestry is shown in Fig. 3. Moreover, rs150717769 is mapped onto the intron in the TRAP1 gene (Table 2). For the sex stratified meta-analysis we found four variants that were significant associated (p < 5e-08) with body shape in women and one in men. These were rs111783937 (TRAP1), rs17867127 (FGF12), rs17867127 (GRM8) , rs75156321 (FGF12), rs17126580 (RNU4ATAC) in females and rs7089940 (TLX1NB) in men. These variants and their related genes were different from those identified for the individual traits. (Fig. S4). We explored other traits these SNPs had been associated with in the Europeans using a PHEWAS approach and the TRAP1 variant was associated with anthropometric trails. This variant based on the Diabetes Epigenome atlas was is involved in regulatory activities as an enahancer in the adipose J o u r n a l P r e -p r o o f tissues. The other variants were inovled in the regulation of the central nervous system as enhancers and bindig sitesin the mammary glands. (Table 2).

Discussion

In this study, we performed a GWASs of PC derived body shape in two different African populations, as well as a meta-analysis of the most similar PCs, which presented several novel genetic loci associated with body shape in individuals of African ancestry. BMI, weight, WC and HC contributed the most in driving body shape patterns and they explained most of the variances for all the six anthropometric traits measured. We identified HNRNPC and GBE1 genetic loci to be associated with body shape in Ugandans and Zulus, respectively. Overall, our meta-analysis of the Ugandan and Zulu GWASs identified TRAP1 as a mutual, additional genetic locus associated with body shape, using the average loadings of the first PC for our six anthropometric traits.

Using average PC from 20 studies, Reid et al. (2016) reported body shape characterised by high loadings for all traits (weight, BMI, WC, WHR, and HC), except height in individuals of European ancestry [START_REF] Ried | A principal component meta-analysis on multiple anthropometric traits identifies novel loci for body shape[END_REF]. Subsequently, in our analysis of continental Africans, we observed that the PC1 had high loadings for weight, BMI, WC and HC, and not in height nor WHR, thus, overall adiposity/body fatness in body shape in Africans is mostly driven by the former four traits. The loadings are in the same direction; meaning that the average PC1 captures inter-individual variation in either increased or decreased BMI, weight, WHR, HC and WC. Interestingly, we observed that the first PC loadings in both Ugandans and Zulus were consistent in weight, BMI, WC and HC, even though the exact variances were different. This observation suggests that the body shape as an index of adiposity has predominantly common measures and features with variability in individuals of African ancestry residing in southern and eastern parts of Africa, as well as Europeans.

On the other hand, the notable differences between Reid et al (2016), and our metaanalysis in the average PC1 are the low loadings of WHR in individuals of African ancestry.

J o u r n a l P r e -p r o o f

Since the majority of our study participants were women, hence, the observed differences that are probably true for women of African ancestry, as they often tend to have wider hips, a lower WHR, and a higher prevalence of steatopygia, than women of European ancestry [START_REF] Danielson | Ethnic variability in bone geometry as assessed by hip structure analysis: findings from the hip strength across the menopausal transition study[END_REF][START_REF] Rush | BMI, fat and muscle differences in urban women of five ethnicities from two countries[END_REF]. Furthermore, we also observed that both the WHR and height in both Ugandans and South African Zulus had low loadings in opposite directions, suggesting the average PC1 does not capture inter-individual variation in either increased or decreased WHR/height. However, the cause is not clear. For example, our second PCs in Ugandans and South African Zulus were driven more by WHR (and then height) and height (and then WHR, both in the opposite direction), respectively, compared to other anthropometric traits, Thus, our individual cohort loadings for our second to sixth PCs were not consistent between Ugandans and South African Zulus to use for meta-analyses.

After linear regression, with blood pressure and blood lipids traits in individuals of African ancestry in Uganda, our first PC was positively associated with adiposity derived body shape. However, in the Zulus, our PC1 was positively associated with blood pressure traits and only LDL; whereas, HDL exerted a protective effect. The observed differences in cardiometabolic trait prediction using PC1 between Ugandans and Zulus may be due to genetic and some environmental exposure differences between these cohorts [START_REF] Yu | Larger Genetic Differences Within Africans Than Between Africans and Eurasians[END_REF][START_REF] Lucchese | Implications of African genetic diversity[END_REF][START_REF] Gomez | Genetic Variation and Adaptation in Africa: Implications for Human Evolution and Disease[END_REF] and future studies are needed to examine this association further.

.Using PC GWAS meta-analysis we identified rs150717769 (p = 4.01e-09) to be associated with body shape in individuals of African ancestry. This variant is mapped in the intergenic region between TRAP1 (Table 2). TRAP1 acts as a key regulator of mitochondrial homeostasis, bioenergetics and metabolism, thus contributing to the regulation of the body shape. Moreover, TRAP1 have been previously identified to be associated with BMI in individuals of diverse ancestry population [START_REF] Turcot | Protein-altering variants associated with body mass index implicate pathways that control energy intake and expenditure underpinning obesity[END_REF]; and TRAP1 was also associated with height in individuals of European ancestry [START_REF] Kichaev | Leveraging Polygenic Functional Enrichment to Improve GWAS Power[END_REF]. Using the GTEx Portal, we observed that TRAP1 is strongly expressed in visceral and subcutaneous adipose tissue, supporting our findings that the average PC1 captures largely adiposity in both Ugandans and South African Zulu. The other variants associated with sex stratified body shape located in GRM8, TLX1NB, FGF12 J o u r n a l P r e -p r o o f were noted to be involved in the regulatory activities in the central nervous system as enhancers. This is similar to what has been reported in related anthropometric GWAS where variants with functinon in the central nervous system have been implicated with obesity.

The main strength of this study is that our PC1 which explains much of anthropometric traits variation was associated with cardiometabolic traits including blood pressure and blood lipids traits in individuals of African ancestry. However, our study was limited by a few SNPs reaching genome-wide significance in both cohorts and our analyses were limited to PC1.

Nevertheless, the other PCs have loadings similar to the single-trait GWASs and they do not explain much of the variation for anthropometric traits.

In summary, we identified five genomic loci to be associated with body shape in individuals of African ancestry that are different from those identified for individual traits in this population group. Our findings suggest that the body shape assessed by PCs provide more detailed information about an individual's body shape, size and adiposity that is not fully captured by individual anthropometric traits. If applied to other highly correlated phenotypes, PCA might reveal novel genetic loci and pathways that have not been identified in single-trait GWAS.

J o u r n a l P r e -p r o o f

study participants in the GPC Ugandan and Zulu cohorts had anthropometric measurements and provided written informed consent to participate in the studies.

Anthropometric traits measurements

At recruitment anthropometric traits including height, weight, WC and HC were measured from all the study participants. Height and weight were measured using the Leicester stadiometer and the Seca 761 class III mechanical flat scales, respectively [START_REF] Gurdasani | Uganda Genome Resource Enables Insights into Population History and Genomic Discovery in Africa[END_REF]. WC and HCs were measured using the non-stretch Seca 201 ergonomic circumference measuring tape. WC was measured at the mid-point between the lower costal margin and the level of the anterior superior iliac crests. HC was measured at the greater trochanter of the femur as previously described [START_REF] Gurdasani | Uganda Genome Resource Enables Insights into Population History and Genomic Discovery in Africa[END_REF]. BMI values were derived from weight over height whilst WHR was derived from WC values over HC values.

Genotype Data

Information on sample collection, genotyping, imputation and quality control procedures were described in detail by previous studies [START_REF] Fatumo | Complimentary Methods for Multivariate Genome-Wide Association Study Identify New Susceptibility Genes for Blood Cell Traits[END_REF][START_REF] Gurdasani | Uganda Genome Resource Enables Insights into Population History and Genomic Discovery in Africa[END_REF] Participant samples from the Uganda GPC were genotyped on the Illumina Human Omni 2.5M Bead Chip array at the Wellcome Trust Sanger Institute (WTSI). The samples were chosen as a subset of the survey population with the most complete phenotype data on the traits measured. Samples from the South African Zulu study were genotyped on the consortium-driven Illumina HumanOmni Multi-Ethnic GWAS/Exome Array (MEGA pre-commercial v1) using the Infinium Assay. Quality control was carried out collectively, with sample QC including filtering for called proportion (< 97%), heterozygosity (> 4SD from mean), sex check fails (F statistic < 0.8 for men, and > 0.2 for women). Sample QC was followed by SNP QC, including filtering for called proportion (< 97%), Hardy Weinberg disequilibrium (p < 1e-06), and relatedness (IBD > 0.90). Samples and variants that did not pass the quality thresholds for the SNP and samples quality control were excluded. Principal Component Analysis (PCA) was carried out to correct for population/ancestry outliers.

Following quality control, phasing and imputation for both cohorts were carried out with SHAPEIT2 using default parameters and IMPUTE2 respectively. Transformation of traits was J o u r n a l P r e -p r o o f carried out uniformly for each cohort to make effect sizes comparable across cohorts, allowing meta-analyses of summary results.

Estimation of PCs loadings and association testing

To capture body shape in its multidimensional structure, we carried out PCA on six anthropometric traits. For each study, we performed a PCA on the standardized residuals of the anthropometric traits and adjusted for age and sex for the Uganda study. Diabetes status in addition to age and sex were corrected for in the South African Zulu study in the combined analysis. In the sex stratified analysis age and diabetes status were corrected for. The standardized residuals were decomposed into six principal components (PCs) according to eigenvectors, and principal loadings. The six resulting PCs is an aggregate and an orthogonal linear representation of individual PCs of the six anthropometric traits. In our analysis loadings is used to describe the weight of each PCs. The loadings and explained variances were comparable between Uganda's and Zulu from South Africa. The association between genetic variants and body shape was tested using the linear mixed model implemented in genomewide efficient mixed-model analysis (GEMMA) [23]. In our analysis, we used PC1 as the phenotype of interest in both Uganda's and Zulu. Furthermore, we adjusted for cryptic relatedness by including a kinship matrix in the model. The kinship matrix was estimated using genetic variants with MAF >1% and pruned to an r2 threshold of 0.5

Meta-Analysis and annotation of genetic variants

Meta-analyses of 8769 participants was performed on the summary statistics of the PC1 of the Ugandan GPC (N = 6178) and South African Zulu cohort (N = 2591), , using the inverse variance weighted fixed-effects model as implemented in Metal. We also lated conducted meta-analysis in men and women for the PC1 bodyshape.This is a free, efficient, and open-source meta-analysis software tool for GWASs analysis [START_REF] Willer | METAL: fast and efficient meta-analysis of genomewide association scans[END_REF]. We used FUMA [START_REF] Watanabe | Functional mapping and annotation of genetic associations with FUMA[END_REF], to annotate, prioritize, visualize, and interpret GWASs results for body shape. Genetic variants with p < 5 x 10 -8 were selected as signals with genome-wide significance.

The SNP2GENE function within FUMA takes GWAS summary statistics as an input and J o u r n a l P r e -p r o o f provides extensive functional annotation in genomic areas identified by lead SNPs, while the GENE2FUNC annotates genes in a biological context [START_REF] Watanabe | Functional mapping and annotation of genetic associations with FUMA[END_REF].

Body shape and cardiometabolic traits

Based on the Kaiser rule, we chose to retain the first and second PCs after PCA. However, in our analysis, we focused on the first PC in Ugandans as it was similar to the first PC in the South African Zulu dataset. We further evaluated whether our PC1 in both cohorts are associated with cardiometabolic traits, we then performed a linear regression analysis in both Ugandans and South African Zulu. In both cohorts, we adjusted for age and sex. We further adjusted for T2D in the South African Zulu dataset. In this analysis, our cardiometabolic traits were blood pressure (systolic and diastolic) and blood lipids (low-density lipoprotein cholesterol, triglycerides, total cholesterol, and high-density lipoprotein cholesterol) traits. The linear regression analysis was performed using the generalised linear model (GLM) and all statistical analyses were performed using the R language [START_REF] Team | A language and environment for statistical computing[END_REF].  Those genetic loci are also known to be associated with Body Mass Index (BMI), height and Waist Circumference.

J o u r n a l P r e -p r o o f

 Our study shows that a standard deviation change in body shape was associated with an increase in blood pressure and blood lipids.

J o u r n a l P r e -p r o o f

Fig 1 .

 1 Fig 1. First PC was used to capture body shape in its multidimensional structure showing loadings for Ugandan and Zulu cohorts. (a) Loading for all participants in the Ugandan cohort, (b) loadings for male-Ugandan cohort (c) loading for female Ugandans. (d) Loading for all participants in the South African Zulu cohort, (e) loadings for male Zulus (f) loadings for female Zulus J o u r n a l P r e -p r o o f

Fig 2 .

 2 Fig 2.PC loadings for Zulu in South Africa using six anthropometric traits to capture body shape in its multidimensional structure. (A) The first PC shows high loading for weight, BMI, WC and HC, (B) The second PC shows high loading for height, (C) The third PC shows high loading for WHR and height, (D) The fourth PC shows high loading WC, but opposite loading for HC, (E) The fifth PC shows high loading for BMI and opposite loadings for WC and HC (F)The sixth PC shows high loading for weight but opposite loading for BMI.J o u r n a l P r e -p r o o f

Fig 3 .

 3 Fig 3. Meta-analysis association between genetic variants and body shape in individuals of African ancestry. (A) Manhattan plot, (B) QQ plot. J o u r n a l P r e -p r o o f

Fig. 4

 4 Fig.4 Locus zoom plot for meta-analysis association between genetic variants and body shape showing rs150717769 mapped within the intergenic region between TRAP1 and DNASE1 J o u r n a l P r e -p r o o f

Fig 2 .

 2 Fig 2.PC loadings for Zulu in South Africa using six anthropometric traits to capture body shape in its multidimensional structure. (A) The first PC shows high loading for weight, BMI, WC and HC, (B) The second PC shows high loading for height, (C) The third PC shows high loading for WHR and height, (D) The fourth PC shows high loading WC, but opposite loading for HC, (E) The fifth PC shows high loading for BMI and opposite loadings for WC and HC (F)The sixth PC shows high loading for weight but opposite loading for BMI.J o u r n a l P r e -p r o o f

Table 1 .

 1 Participant characteristics

	Variables	Ugandan study (n = 6,407)	South African Zulu (n = 2,595)
		median (IQR)	median (IQR)
	Female, n (%)	3,660 (57.1)	1,919 (73.8)
	Age, years	29.8 (17.8-46.1)	52.0 (39.0-60.0)
	BMI, kg/m 2	20.7 (18.8-23.0)	31.4 (26.1-36.9)
	Height, m 2	1.57 (1.51-1.63)	1.60 (1.55-1.65)
	Weight, kg	52.0 (45.0-60.0)	80.0 (67.5-93.6)
	WHR	0.84 (0.80-0.88)	0.92 (0.85-1.00)
	HC, cm	89.0 (83.5-94.5)	109 (100-120)
	WC, cm	74.0 (69.1-79.5)	103 (91.0-113)
	IQR; interquartile range, BMI; body mass index, WHR; waist-hip-ratio, HC; hip circumference, WC; weight
	circumference

J o u r n a l P r e -p r o o f

Acknowledgements

SF is an international Intermediate fellow funded by the Wellcome Trust grant (220740/Z/20/Z) at the MRC/UVRI and LSHTM. TC is an international training fellow supported by the Wellcome Trust grant (214205/Z/18/Z). MN is a Makerere University Non-Communicable Diseases (MakNCD) Fellow. MN acknowledges the support of the Makerere University Non-Communicable Diseases (MakNCD) Research Training Program supported by the Fogarty International Center (FIC) of the National Institutes of Health (NIH) under award number 1D43TW011401-01. MN and SF received previous support from the Makerere University-Uganda Virus Research Institute Centre of Excellence for Infection and Immunity Research and Training (MUII). MUII is supported through The Developing Excellence in Leadership, Training and Science (DELTAS) Africa Initiative (grant 107743). The DELTAS Africa Initiative is an independent funding scheme of the African Academy of Sciences (AAS), Alliance for Accelerating Excellence in Science in Africa (AESA), and supported by the New Partnership for Africa's Development Planning and Coordinating Agency (NEPAD Agency) with funding from the Wellcome Trust (107743) and the UK government. TM is a RHDGen PhD fellow supported by the Wellcome Trust, the University of Cape Town, and the inaugural Bongani Mayosi UCT-PHRI Scholarship (McMaster University).

Methods

Study Population

The Medical Research Council (MRC) United Kingdom, together with the Uganda Virus Research Institute (UVRI) launched the General Population Cohort in the 1980s with the primary aim of investigating the epidemiology of Human Immunodeficiency Virus (HIV) infection in Uganda's. In 2011, the University of Cambridge, Wellcome Sanger Institute (WSI), in collaboration with the MRC/UVRI began the GPC round 22 study to investigate the genetics and epidemiology of communicable and non-communicable diseases in Uganda's, using both populations genetic and epidemiological approaches. A total of 4,778 Uganda's from the GPC round 22 study and 1,629 individuals from the Uganda 2,000 Genomes project (UG2G) [START_REF] Gurdasani | Uganda Genome Resource Enables Insights into Population History and Genomic Discovery in Africa[END_REF][START_REF] Fatumo | Complimentary methods for multivariate genome-wide association study identify new susceptibility genes for blood cell traits[END_REF], making up a total of 6,407 individuals were used as the discovery cohort.

Conflicts of interest

None

J o u r n a l P r e -p r o o f


Paulo C T Souza
email: paulocts@gmail.com

Riccardo Alessandri

Jonathan Barnoud

Sebastian Thallmair

Ignacio Faustino

Fabian Grünewald

Ilias Patmanidis

Haleh Abdizadeh

Bart M H Bruininks

Tsjerk A Wassenaar

Fabian Grunewald

Peter C Kroon

Josef Melcr

Vincent Nieto

Valentina Corradi

Hanif M Khan

Jan Domański

Matti Javanainen

Hector Martinez-Seara

Nathalie Reuter

Robert B Best

Ilpo Vattulainen

Luca Monticelli

Xavier Periole

D Peter Tieleman

Alex H De Vries

Siewert J Marrink

Martini 3: a general purpose force field for coarse-grained molecular dynamics

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The molecular dynamics (MD) simulation technique has become an indispensable tool in natural sciences, offering a spatio-temporal resolution unmatched by any experimental technique [START_REF] Bottaro | Biophysical experiments and biomolecular simulations: A perfect match?[END_REF] . A major bottleneck of MD is the limited time and length scales that are accessible. To overcome this limitation, coarse-grained (CG) models representing groups of atoms by effective beads, have achieved widespread use [START_REF] Ingólfsson | The power of coarse graining in biomolecular simulations[END_REF] . The Martini model is among the most popular CG models in the field of biomolecular simulation, due to its easy-touse building block principle. Martini relies on a four-to-one mapping scheme (i.e., on average four heavy atoms and associated hydrogens are mapped into one CG bead), and has been parametrized using a top-down approach with thermodynamic partitioning data as the main target [START_REF] Marrink | Coarse Grained Model for Semiquantitative Lipid Simulations[END_REF][START_REF] Marrink | The MARTINI force field: coarse grained model for biomolecular simulations[END_REF] . Non-bonded interactions between neutral beads of Martini are solely described by Lennard-Jones (LJ) potentials, while charged beads also include Coulombic interactions. The interaction strength of the LJ potential (i.e. its well depth) is used to discriminate between different levels of polarity of the CG beads. The model features four main classes of CG bead types, denoted C, N, P, and Q representing non-polar, intermediately polar, polar, and charged chemical groups, respectively [START_REF] Marrink | The MARTINI force field: coarse grained model for biomolecular simulations[END_REF] . Sub-labels are used to make a further distinction within a class in terms of degree of polarity or hydrogen donor/acceptor capabilities.

In principle, all beads have the same size, denoted as regular (R) beads. By way of exception, special small (S) beads were introduced to model ring-like compounds for which a four-to-one mapping scheme is inadequate [START_REF] Marrink | The MARTINI force field: coarse grained model for biomolecular simulations[END_REF] . To reproduce correct stacking and hydrogen bonding distances between nucleotides, even smaller tiny (T) beads were found necessary [START_REF] Uusitalo | Coarse-Grained Force Field: Extension to DNA[END_REF] .

Parametrization of the cross-interactions between S-and T-beads with Rbeads, however, was done on an ad-hoc basis.

The Martini force field is used in a wide range of applications in diverse fields including structural biology [START_REF] Abellón-Ruiz | Structural basis for maintenance of bacterial outer membrane lipid asymmetry[END_REF][START_REF] Yen | PtdIns(4,5)P2 stabilizes active states of GPCRs and enhances selectivity of G-protein coupling[END_REF][START_REF] Van Eerden | Exchange pathways of plastoquinone and plastoquinol in the photosystem II complex[END_REF] , biophysics [START_REF] Vögele | Hydrodynamics of Diffusion in Lipid Membrane Simulations[END_REF][START_REF] Agostino | A tethering complex drives the terminal stage of SNARE-dependent membrane fusion[END_REF] , biomedicine [START_REF] Jeena | Mitochondria localization induced self-assembly of peptide amphiphiles for cellular dysfunction[END_REF] , nanotechnology [START_REF] Jiang | Subnanometre ligand-shell asymmetry leads to Janus-like nanoparticle membranes[END_REF][START_REF] Maingi | Stability and dynamics of membrane-spanning DNA nanopores[END_REF] , and materials design [START_REF] Frederix | Exploring the sequence space for (tri-)peptide self-assembly to design and discover new hydrogels[END_REF][START_REF] Bochicchio | Into the dynamics of a supramolecular polymer at submolecular resolution[END_REF] . With its growing use, however, a number of shortcomings of the Martini model have recently been identified.

One of the most important problems is the observation that certain molecules tend to interact too strongly. This has been reported for proteins and carbohydrates in solution, as well as for membrane embedded proteins [START_REF] Stark | Toward optimized potential functions for protein-protein interactions in aqueous solutions: osmotic second virial coefficient calculations using the MARTINI coarse-grained force field[END_REF][START_REF] Javanainen | Excessive aggregation of membrane proteins in the Martini model[END_REF][START_REF] Schmalhorst | Overcoming the Limitations of the MARTINI Force Field in Simulations of Polysaccharides[END_REF] .

The origin lies among others in small but systematic deviations in packing and intermolecular interactions [START_REF] Alessandri | Pitfalls of the Martini Model[END_REF] . Besides, the coverage of chemical space for broader applications was uneven, and in some cases, such as selectivity of nucleobase pairing [START_REF] Uusitalo | Coarse-Grained Force Field: Extension to DNA[END_REF][START_REF] Uusitalo | Martini Coarse-Grained Force Field: Extension to RNA[END_REF]

RESULTS

Re-parametrization of the beads. In Martini 3, the new parametrization strategy was based on the construction of prototype models of polar and nonpolar molecules in all three Martini resolutions. Self-(R-R, S-S, and T-T) and cross-interactions (R-S, R-T, and S-T) of the different bead sizes were optimized to be well-balanced (Supplementary Notes, sections B1 and B2).

In terms of chemical types, the beads were separated into three blocks: organic, ions, and water (Supplementary Notes, section A1). The organic (containing P-, N-, and C-beads) and ion (Q-beads) blocks have been subjected to independent parametrizations, where different trends in self- The improved interaction balance between regular and smaller bead types is illustrated by the close to ideal mixing behavior of pure solvents composed of molecules mapped at different resolutions (Fig. 1A). Integration of radial distribution functions, defined as Kirkwood-Buff integrals (G ij ), are used here to quantify the degree of miscibility of the multi-resolution liquid water model (Fig. 1B). Theoretically, pair differences in Kirkwood-Buff integrals (∆G ij ) should be equal to zero for all i,j pairs in ideal mixtures [START_REF] Ben-Naim | Molecular theory of solutions[END_REF][START_REF] Ploetz | Kirkwood-Buff integrals for ideal solutions[END_REF] while real mixtures that closely approach ideal behavior (like benzene-toluene) show values around ± 1 cm 3 /mol [START_REF] Ploetz | Kirkwood-Buff integrals for ideal solutions[END_REF] Another benefit of the recalibrated interactions is the disappearance of the artificially large desolvation free energy barriers that contribute to slow dissociation processes of the previous Martini 2 models. The problem, that was initially observed in dimerization of nucleobases [START_REF] Uusitalo | Coarse-Grained Force Field: Extension to DNA[END_REF][START_REF] Alessandri | Pitfalls of the Martini Model[END_REF] , is thus solved, as highlighted by the comparison of Martini 2 and Martini 3 PMFs between cytosine and guanine (left panel of Fig. 1D). Note that there is room for further improvement, as the free energy minima of the CG PMF profiles with Martini 3 are shifted relative to the all-atom profiles because the bead sizes representing nitrogen-containing groups are not optimal to reproduce hydrogen bonding distances. In addition, the difference between C-G and G-G base pairs is not as large as in the atomistic case (~20 kJ/mol). However, it is still large enough (~8 kJ/mol) to provide specificity.

The proper balancing of R-, S-, and T-beads in Martini 3 also implies that the mapping choice of an arbitrary molecule to its Martini representation is now better defined. S-and T-beads are not only suited to represent ring-like compounds, but also used for cases involving 3-to-1 and 2-to-1 mapping of linear and branched chemical groups (Supplementary Notes, section C1).

Covering the chemical space with new beads and labels. Together with a thorough revision of the interaction strengths, in Martini 3, we extend the number of chemical bead types and the ability to modify the bead properties depending on the chemical details of the underlying moieties. Each P-, N-, and C-class now has 6 bead types with different degrees of polarity, which enables a more precise definition of different chemical groups by assigning them to certain bead types. Additionally, we introduce a new X-class of beads to model halo-compounds (Supplementary Notes, section A1). In the previous version of Martini, some of the bead types were already subclassified according to their ability to act as hydrogen bond donor, acceptor, or both. This property can be now attributed to all bead types of intermediate or polar nature (N-or P-class). The effective interaction strength between donor and acceptor pairs is increased, whereas donor-donor and acceptor-acceptor pairs are weakened (Supplementary Notes, section A4). For example, Martini 3 correctly reproduces the trends in hydrogen bond-based pairing of nucleobases [START_REF] Uusitalo | Coarse-Grained Force Field: Extension to DNA[END_REF][START_REF] Uusitalo | Martini Coarse-Grained Force Field: Extension to RNA[END_REF] without the use of special-purpose beads specifically for nucleobases (right panel of Fig. 1D). Note that chemical groups that can act as both donor and acceptor at the same time are always represented by the pure beads of the P-and N-class in Martini 3.

Next to the fine-tuning based on hydrogen bonding capabilities, we introduce the possibility to change the interactions based on the electronic polarizability.

Depending on inductive or conjugate effects caused by chemical functionalization, non-polar molecules can be polarized, i.e., they can acquire an electron-donor (or "enriched", label "e") or electron-acceptor (or "vacancy", label "v") character, which can promote preferential interactions. Polarizable groups in Martini 3 can be distinguished through the label "e/v" which can only be applied to the C-and X-class. A nice example of their application is the strong and specific interaction between electron donor and electron acceptor aromatic rings in aedamers, a class of molecules that have been studied extensively in the context of biomimetic folding and self-assembly [START_REF] Zych | Synthesis and Conformational Characterization of Tethered, Self-Complexing 1,5-Dialkoxynaphthalene/1,4,5,8-Naphthalenetetracarboxylic Diimide Systems[END_REF][START_REF] Gabriel | Aromatic oligomers that form hetero duplexes in aqueous solution[END_REF] . The use of "e/v" allows Martini 3 to capture the preferential interaction between 1,5-dialkoxynaphthalene (DAN) and naphthalene diimide (NDI) (left panel of Fig. 2A) experimentally observed via NMR titration [START_REF] Zych | Synthesis and Conformational Characterization of Tethered, Self-Complexing 1,5-Dialkoxynaphthalene/1,4,5,8-Naphthalenetetracarboxylic Diimide Systems[END_REF] and atomistic simulation data. Self-assembly of amide-linked tetramers shows preferential formation of alternating stacks of DAN and NDI, which is also measured by NMR and isothermal titration calorimetry investigations [START_REF] Gabriel | Aromatic oligomers that form hetero duplexes in aqueous solution[END_REF] . On top of hydrogen bonding and electron polarization labels, all beads can have their self-interaction finetuned by other sub-labels (as further described in the Supplementary Notes, section A4).

Chemical groups carrying monovalent charges +1/-1 are represented in

Martini by the class of Q-beads (Supplementary Notes, section A1). The original Martini model only considers monovalent ions, and was solely optimized for regular bead sizes which represented small ions and their first hydration shell. In Martini 3, charged groups can have either R-, S-, or T-size.

The tiny size category allows modeling of small, bare ions, enabling applications that involve ion binding where (part of) the hydration shell is lost.

This feature is exemplified by the binding of sodium ions (represented by a charged tiny bead) to a buried small cavity localized in the core of the adenosine A 2A receptor (Fig. 2B). X-ray crystallographic [START_REF] Liu | Structural basis for allosteric regulation of GPCRs by sodium ions[END_REF] and ligand binding assays [START_REF] Gao | Allosteric modulation of A(2A) adenosine receptors by amiloride analogues and sodium ions[END_REF] confirm the importance of sodium ions for the structure and for the allosteric modulation of the A 2A receptor. Note that an extensive validation of the lipid models in Martini 3 was performed to allow simulations of transmembrane and peripheral membrane proteins (see Supplementary Results, section F1).

In addition to the smaller sizes, the Q-class was also expanded to five bead types, following the classical Hofmeister series trend [START_REF] Okur | Beyond the Hofmeister Series: Ion-Specific Effects on Proteins and Their Biological Functions[END_REF][START_REF] Dupont | Overview of the Effect of Salts on Biphasic Ionic Liquid/Water Solvent Extraction Systems: Anion Exchange, Mutual Solubility, and Thermomorphic Properties[END_REF] Br -(SQ4 -bead) or (partially) PF 6 -(Q1 -bead) from the ionic liquid phase. In the case of the biphasic system using trihexyltetradecylphosphonium bromide

([P 66614 ][Br]
), direct comparison to experimental data shows good agreement for the anion transfer trends [START_REF] Dupont | Overview of the Effect of Salts on Biphasic Ionic Liquid/Water Solvent Extraction Systems: Anion Exchange, Mutual Solubility, and Thermomorphic Properties[END_REF][START_REF] Naert | Ionic liquid ion exchange: Exclusion from strong interactions condemns cations to the most weakly interacting anions and dictates reaction equilibrium[END_REF] . The new Q-bead types also impact biologically relevant systems, as exemplified by preferential cation-π interaction between choline groups (Q1 + bead) of phosphatidylcholine lipids and aromatic residues of the Bacillus thuringiensis phosphatidylinositolspecific phospholipase C (BtPI-PLC). In the previous version, such specific interaction between soft ions and aromatic molecules were solely included in the recently updated polarizable Martini implementation [START_REF] Khan | Capturing Choline-Aromatics Cation-πInteractions in the MARTINI Force Field[END_REF] . However, in Martini (FraC) nanopore inserted in a lipid bilayer (Fig. 3A). The pore of FraC remains open over the whole trajectory in Martini 3, as indicated by X-ray crystallography [START_REF] Tanaka | Structural basis for self-assembly of a cytolytic pore lined by protein and lipid[END_REF] and electro-osmotic flow assays 32 , while in Martini 2 it is closed.

Another example of accurate packing is the stacking predictions of thiophene derivatives in bulk heterojunction solar cells composed of poly(3-hexylthiophene) (P3HT) and phenyl-C61-butyric acid methyl ester (PCBM) (Fig. 3B). The morphology of these organic solar cells is a determinant for highefficiency devices [START_REF] Alessandri | Bulk Heterojunction Morphologies with Atomistic Resolution from Coarse-Grain Solvent Evaporation Simulations[END_REF] . The scattering profiles computed with Martini 3 show improved agreement with Martini 3 in relation to P3HT lamellar (peak around q≈0.45 Å -1 ) and stacking (q≈1.65 Å -1 ) experimental distances [START_REF] Alessandri | Bulk Heterojunction Morphologies with Atomistic Resolution from Coarse-Grain Solvent Evaporation Simulations[END_REF][START_REF] Chiu | Simultaneous use of small-and wide-angle X-ray techniques to analyze nanometerscale phase separation in polymer heterojunction solar cells[END_REF] .

The use of bonds based on molecular volume and the appropriate choice of chemical bead types, sizes, and labels also controls the interaction density of the model, which has an important impact on the strength of collective interactions between molecules [START_REF] Alessandri | Pitfalls of the Martini Model[END_REF] . In order to test to what extent the changes in non-bonded and bonded interactions reduce the over-estimated aggregation of proteins, we performed extensive simulations comprised of solutions of soluble proteins as well as membrane embedded proteins. These systems were simulated under conditions in which proteins do not aggregate and, preferentially, stay as monomers. For soluble proteins (Fig. 3C),

qualitative tests were performed with the headpiece domain of chicken villin [START_REF] Petrov | Are Current Atomistic Force Fields Accurate Enough to Study Proteins in Crowded Environments?[END_REF] , and the modified and mutated cellulose-binding domain from Cellulomonas fimi (EXG-CBM), which is an example of a protein completely free of charged side chains that can maintain solubility, stability, and function [START_REF] Højgaard | A Soluble, Folded Protein without Charged Amino Acid Residues[END_REF] . Trends are well-captured in Martini 3, with both proteins mainly staying as monomers in pure water (with only counter-ions to neutralize the system in the case of villin). The villin headpiece showed salting-in behavior (i.e. less aggregation) under addition of 0.4M of NaCl, which was also observed for certain soluble proteins at low ionic strengths [START_REF] Ruckenstein | Effect of salts and organic additives on the solubility of proteins in aqueous solutions[END_REF] . On the other hand, EXG-CBM only showed salting-out behavior (i.e. more aggregation), which was expected based on experimental data [START_REF] Højgaard | A Soluble, Folded Protein without Charged Amino Acid Residues[END_REF] . In contrast, both proteins aggregate in Martini 2, forming a single and stable aggregate during the simulation.

Polyleucine (K 2 -L 26 -K 2 ) was selected to evaluate the aggregation propensity in membranes. Experimental evidence with this transmembrane (TM) protein model indicates a preference for the monomeric state in a bilayer environment [START_REF] Zhou | Interhelical hydrogen bonding drives strong interactions in membrane proteins[END_REF][START_REF] Zhou | Polar residues drive association of polyleucine transmembrane helices[END_REF][START_REF] Grau | The role of hydrophobic matching on transmembrane helix packing in cells[END_REF] . Both Martini 2 and 3 show that the hydrophobic mismatch between TM length and membrane thickness can play a role in the aggregation, with Martini 3 showing a higher percentage of the monomeric state (Fig. 3D). To quantitatively evaluate the strength of protein-protein interactions in a membrane environment, we also considered the dimerization of four selected transmembrane (TM) helices: the TM domains of the receptor tyrosine kinases EphA1 and ErbB1; the red blood cell protein glycophorin A (GpA); as well as the well-known model peptide WALP23 (left panel of Fig. 3E). For EphA1 and ErbB1, experimental dimerization free energies in a membrane environment have been estimated using Förster resonance energy transfer (FRET) [START_REF] Chen | Energetics of ErbB1 transmembrane domain dimerization in lipid bilayers[END_REF][START_REF] Artemenko | Transmembrane domain of EphA1 receptor forms dimers in membrane-like environment[END_REF] . For GpA, dimerization free energies range from around -15 kJ/mol (in various cell membrane environments) [START_REF] Sarabipour | Glycophorin A transmembrane domain dimerization in plasma membrane vesicles derived from CHO, HEK 293T, and A431 cells[END_REF][START_REF] Chen | Measuring the Energetics of Membrane Protein Dimerization in Mammalian Membranes[END_REF] to -31.5 kJ/mol (GALLEX assay in E. coli inner membranes) [START_REF] Nash | De novo design of transmembrane helix-helix interactions and measurement of stability in a biological membrane[END_REF][START_REF] Finger | The Stability of Transmembrane Helix Interactions Measured in a Biological Membrane[END_REF] or -50.6 kJ/mol (steric trap in POPC bilayers) [START_REF] Hong | Method to measure strong protein-protein interactions in lipid bilayers using a steric trap[END_REF] . WALP peptides have been characterized thoroughly during the past two decades, including their self-association [START_REF] Sparr | Self-association of transmembrane α-helices in model membranes: Importance of helix orientation and role of hydrophobic mismatch[END_REF] . For each one of the four peptide dimers, we compared experimental dimerization free energies with the free energies predicted by the Martini 2 and Martini 3 model predictions. Martini 3 shows not only to be able to capture the correct trends, but also to quantitatively reproduce the experimental affinities. The binding mode also becomes improved as highlighted for GpA (right panel of Fig. 3E).

The GpA homodimer structure with Martini 3 closely resembles experimental results obtained with NMR spectroscopy and crystallography [START_REF] Mackenzie | Transmembrane helix dimer: Structure and implications[END_REF][START_REF] Trenker | Crystal Structure of the Glycophorin A Transmembrane Dimer in Lipidic Cubic Phase[END_REF][START_REF] Domański | Balancing Force Field Protein-Lipid Interactions To Capture Transmembrane Helix-Helix Association[END_REF] .

In summary, for both soluble and transmembrane proteins, we observed that the Martini 3 models are in much better agreement with experimental data than before. Another advantage of the current Martini 3 protein model is that the default bead type representing the protein backbone in Martini 3 (a regular P2 bead) no longer depends on the secondary structure. In addition, the representation of protein flexibility can now be improved by the use of Gō-like models [START_REF] Souza | An Allosteric Pathway in Copper, Zinc Superoxide Dismutase Unravels the Molecular Mechanism of the G93A Amyotrophic Lateral Sclerosis-Linked Mutation[END_REF] .

DISCUSSION

In this paper we have described the new version of the Martini force field, which shows numerous improvements in relation to the previous version.

However, inherent limitations to the process of coarse-graining, related to transferability and representability problems [START_REF] Brini | Systematic coarse-graining methods for soft matter simulations-a review[END_REF][START_REF] Foley | The impact of resolution upon entropy and information in coarse-grained models[END_REF][START_REF] Wagner | On the representability problem and the physical meaning of coarse-grained models[END_REF] are still part of the model. An important drawback is the limited structural detail, that is a consequence of representing multiple atoms with isotropic interaction sites. This is most noticeable for the Martini water model, which represents four water molecules with a single LJ site and will certainly not capture any of the higher order structural correlations of real water. The role of explicit water in a CG model such as Martini is mostly to provide a good solvent for polar compounds resulting in realistic partitioning. For applications requiring finer details, structure-based CG models are more suitable [START_REF] Wörner | Direct route to reproducing pair distribution functions with coarse-grained models via transformed atomistic cross correlations[END_REF][START_REF] Noid | Multiscale coarsegraining and structural correlations: Connections to liquid-state theory[END_REF] . Another fundamental limitation is the entropy-enthalpy compensation. The loss of internal degrees of freedom for groups of atoms represented by a CG bead inevitably reduces the entropy of the system. Since the Martini force field is based on reproducing free energies, this requires a concomitant reduction in the enthalpy. As consequence, inaccurate entropy-enthalpy balance affects the temperature dependence of several properties and reduces the transferability to different state points. To probe transferability, we performed temperature dependent calculations for a number of solvent systems as well as lipid membranes (Supplementary Results F). Temperature dependent properties, like the heat expansion coefficient and heat capacity of water and n-octane, are very well captured by Martini 3, but this is not true for the hydrophobic effect, that shows the opposite trend with respect to atomistic models, in line with previous findings [START_REF] Wu | Driving force for the association of hydrophobic peptides: The importance of electrostatic interactions in coarsegrained water models[END_REF] . Note that the use of higher-resolution S-or Tparticles does not remedy this problem, as these bead types were parameterized mainly to be compatible with the regular (R-type) beads and should be used primarily to represent parts of the system that cannot be adequately mapped with R-particles. Potential improvements with respect to the temperature transferability of our model could make use of environment dependent potentials [START_REF] Jin | Temperature and Phase Transferable Bottomup Coarse-Grained Models[END_REF] or CG beads with embedded sites, such as the polarizable water models [START_REF] Yesylevskyy | Polarizable water model for the coarse-grained MARTINI force field[END_REF][START_REF] Michalowsky | A refined polarizable water model for the coarse-grained MARTINI force field with longrange electrostatic interactions[END_REF] , where incorporation of quadrupole moment might be required [START_REF] Wu | Driving force for the association of hydrophobic peptides: The importance of electrostatic interactions in coarsegrained water models[END_REF] . Bottom-up CG models that are derived with minimization of the information loss [START_REF] Foley | The impact of resolution upon entropy and information in coarse-grained models[END_REF] as parameterization target might also perform better.

For a more in depth discussion of these and related issues with respect to the Martini coarse-graining philosophy, we refer to previous papers [START_REF] Marrink | Perspective on the Martini model[END_REF][START_REF] Bruininks | A Practical View of the Martini Force Field[END_REF] .

Keeping these limitations in mind, Martini 3 offers a versatile and easy to use generic force field to simulate molecular processes at a semi-quantitative level of accuracy. In relation to the previous model, the excessive over-estimated aggregation [START_REF] Alessandri | Pitfalls of the Martini Model[END_REF] is substantially reduced. We expect that Martini 3 will allow for more realistic predictions of protein interactions, as well as more accurate simulations of molecular systems in general. The increased number of bead types and interaction levels makes the model even more versatile, covering a larger part of chemical space with appropriate building blocks. Based on this new foundation, further optimizations for different classes of molecules are currently ongoing, including the use of Gō-potentials to alleviate limitations of protein conformational flexibility, a re-optimization of the bonded potentials of lipids and other biomolecular classes, as well as a compatible polarizable water model for applications requiring more realistic screening of electrostatic interactions. Finally, we foresee new application horizons for the Martini model in the field of materials science [START_REF] Liu | Enhancing Molecular n-Type Doping of Donor-Acceptor Copolymers by Tailoring Side Chains[END_REF][START_REF] Vazquez-Salazar | Martini coarse-grained models of imidazolium-based ionic liquids: from nanostructural organization to liquid-liquid extraction[END_REF] and high-throughput drug design [START_REF] Souza | Protein-ligand binding with the coarse-grained Martini model[END_REF] .

ONLINE CONTENT

Methods, including statements of data and code availability, supplementary information, and any other associated content and references, are available in the online version of the paper at https://doi.org/10.1038/sXXXXX-XXX-XXXX-X. bootstrapping while experimental data was obtained in the literature [START_REF] Chen | Energetics of ErbB1 transmembrane domain dimerization in lipid bilayers[END_REF][START_REF] Artemenko | Transmembrane domain of EphA1 receptor forms dimers in membrane-like environment[END_REF][START_REF] Sarabipour | Glycophorin A transmembrane domain dimerization in plasma membrane vesicles derived from CHO, HEK 293T, and A431 cells[END_REF][START_REF] Chen | Measuring the Energetics of Membrane Protein Dimerization in Mammalian Membranes[END_REF][START_REF] Nash | De novo design of transmembrane helix-helix interactions and measurement of stability in a biological membrane[END_REF][START_REF] Finger | The Stability of Transmembrane Helix Interactions Measured in a Biological Membrane[END_REF][START_REF] Hong | Method to measure strong protein-protein interactions in lipid bilayers using a steric trap[END_REF][START_REF] Sparr | Self-association of transmembrane α-helices in model membranes: Importance of helix orientation and role of hydrophobic mismatch[END_REF] . In the case of GpA, error was estimated by the mean absolute error of four independent experimental data 43-47 . A comparison between experimental and simulated binding modes of GpA is highlighted in the right panel. The experimental structure was taken from solution NMR in micelles (PDB code AFO) [START_REF] Mackenzie | Transmembrane helix dimer: Structure and implications[END_REF] . Near identical experimental structures were obtained by ssNMR in nanodiscs and Xray crystallography in a lipid cubic phase [START_REF] Mackenzie | Transmembrane helix dimer: Structure and implications[END_REF][START_REF] Trenker | Crystal Structure of the Glycophorin A Transmembrane Dimer in Lipidic Cubic Phase[END_REF][START_REF] Domański | Balancing Force Field Protein-Lipid Interactions To Capture Transmembrane Helix-Helix Association[END_REF] .

mapping was inspired by the previous Martini model [START_REF] Marrink | The MARTINI force field: coarse grained model for biomolecular simulations[END_REF][START_REF] Wassenaar | Computational Lipidomics with insane : A Versatile Tool for Generating Custom Membranes for Molecular Simulations[END_REF][START_REF] Melo | Parameters for Martini sterols and hopanoids based on a virtual-site description[END_REF][START_REF] López | Force Field Parameters for Glycolipids[END_REF] , but now following the Martini 3 rules for mapping and also with adaptations in the bonded parameters inspired by the "extensible model" of Carpenter et al. [START_REF] Carpenter | Capturing Phase Behavior of Ternary Lipid Mixtures with a Refined Martini Coarse-Grained Force Field[END_REF] .

2) General setup for CG MD simulations and Analysis

Settings for the CG simulations followed, in general, the "new" Martini set of simulation parameters [START_REF] De Jong | Martini straight: Boosting performance using a shorter cutoff and GPUs[END_REF] using the leap-frog algorithm [START_REF] Hockney | Quiet high-resolution computer models of a plasma[END_REF] for integrating the equations of motion. The Verlet neighbor search algorithm [START_REF] Páll | A flexible algorithm for calculating pair interactions on SIMD architectures[END_REF] is used to update the neighbor list every 20 steps with a buffer tolerance of 0.005 kJ•mol -1 •ps -1 .

For the Lennard-Jones terms, we used a cutoff scheme with a value of 1.1 nm and the Verlet cutoff scheme [START_REF] Verlet | Computer 'experiments' on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF] for the potential-shift. Long range electrostatic interactions were treated with reaction field [START_REF] Tironi | A generalized reaction field method for molecular dynamics simulations[END_REF] or Particle Mesh Ewald (PME) [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF] , with relative permittivity set to =15 and a cutoff value of 1.1 nm. Reaction field was used for most of the systems, except the ones explicitly mentioning PME. Periodic boundary conditions were used in all three dimensions. For the production simulations, the velocity rescaling thermostat [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF]

3) Parameter calibration, tests and validation

In order to parametrize the LJ parameters of single beads and also test the Martini 3 CG models, many molecular systems and methods were used in this work. The overall iterative approach was not based in rigorous separation of calibration and validation groups. As Martini is based on pair interactions, it is hard to find simple systems that cover enough points in the interaction matrix for all bead size combinations. So, complex systems are not only used for validation but can be part of the calibration. The tests performed were separated in "tiers", which represent systems with different level of complexity.

In "tier 0", isolate beads and simple-molecules are mainly used for calibration of LJ parameters, with balance of different bead sizes and thermodynamics data (e.g. liquid-liquid partitioning and miscibility) used as main targets. In the intermediate "tier 1", bilayer properties are checked, together with qualitative tests, applied to systems like soluble and transmembrane proteins. These qualitative tests are designed as "yes-or-no" questions to evaluate the overall quality of the force field. At the same time, some points in the interaction matrix were also tested and fine-tuned here. In the final "tier 2", quantitative

3 ,

 3 the different Q-bead types allow easier (implicit) incorporation of such interactions without the need for additional partial charges.On top of the new chemical types, all Q-beads can use the hydrogen-bonding labels (called in this case "p/n"). They represent organic charged molecules or fragments, such as R-CH 2 -COO -and R-CH 2 -NH 3 + , and also introduce modifications in the Hofmeister trends of the pure Q-beads (Supplementary Notes, section A4). Positively charged hydrogen donors ("p" label) interact more strongly with non-polar beads, as expected in cation-π interactions. On the other hand, negatively charged hydrogen acceptors ("n" label), have stronger interactions with neutral polar beads, mimicking the stronger hydrogen-bonds with anions. To complete the ion block, we explicitly include a new D-bead for divalent ions (such as Ca 2+ ), which are typically hard ions.Improving packing and protein-protein interactions. Another change inphilosophy with respect to the previous Martini models is the definition of bonded interactions. Instead of using the center of mass of the mapped chemical groups to define the geometry of the molecule, we now use a sizeshape concept aimed at preserving the volume of molecules in comparison to all-atom reference structures. This choice and the proper use of Martini 3 bead sizes (Supplementary Notes, sections C1 and C2) lead to more realistic molecular packing. As a consequence, the hydration of protein pockets and channels is improved, as illustrated by the Fragaceatoxin C

Figure 1 :

 1 Figure 1: Rebalancing R-, S-, and T-beads -(A) Snapshots of simulation boxes containing mixtures of dodecane and water in three resolutions. (B) Radial distribution functions (g ij ) for all bead combinations in the multi-resolution mixture of water. (C) Water-oil transfer free energies (∆G) computed for around 260 data points using Martini 3. (D) Hydrogen bonding potential of mean force (PMF) between nucleobases. On the left, comparison between Martini 2 and 3 for the cytosine-guanine base pair. On the right, comparison of the cytosine-guanine (C-G) and guanine-guanine (G-G) base pairs using Martini 3. In both plots, CHARMM andAMBER atomistic data are also reported[START_REF] Uusitalo | Coarse-Grained Force Field: Extension to DNA[END_REF] for comparison. Errors are estimated with bootstrapping and displayed as transparent shades. In the case of Martini, errors are smaller than 0.1 kJ/mol, and hence are not visible in the graphs.

Figure 2 :

 2 Figure 2: New chemical bead types, sub-labels, and applications -(A) Self-assembly of aedamers. The left panel shows the dimerization free energies (∆G dim ) of pegylated monomers of 1,5-dialkoxynaphthalene (DAN) and naphthalene diimide (NDI). Errors are estimated with bootstrapping. The right panel shows the self-assembled duplex dimer formed by amide-linked tetramers of NDI (green) and DAN (orange). (B) As indicatedby X-ray crystallography 25 , sodium ions (charged TQ5 bead) can bind to a buried small cavity in the core of the adenosine A 2A receptor. (C) Charged Q-beads in Martini 3 follow the classical Hofmeister series, as exemplified by the anion transfer between salt aqueous solutions and organophosphonium-based ionic liquids (right panel). Errors in the average anion transfer percentage are estimated by block averaging. (D) Preferential cation-π interaction between choline groups (Q1 bead) of phosphatidylcholine lipids and aromatic residues of the Bacillus thuringiensis phosphatidylinositol-specific phospholipase C (BtPI-PLC). The depth of insertion of each amino acid of BtPI-PLC is in very good agreement with the insertion obtained from an atomistic MD simulation 30 .

Figure 3 :

 3 Figure 3: Improving packing, cavities and reducing protein stickiness -(A) Hydration of Fragaceatoxin C (FraC) nanopore inserted in a lipid bilayer. (B) Scattering profiles and a Martini 3 snapshot of a bulk heterojunction morphology of poly(3-hexyl-thiophene) (P3HT, in red) and phenyl-C61-butyric acid methyl ester (PCBM, in blue) formed after solvent evaporation and annealing simulations. I(q) corresponds to scattering intensity and q is the reciprocal space vector. (C) Aggregation levels of the soluble proteins villin headpiece and the modified EXG-CBM in different salt concentrations. (D) Aggregation levels of polyleucine helices in POPC and DLPC bilayers. Errors in the average monomer percentage of (C) and (D) are estimated by block averaging. (E) Dimerization of transmembrane helices. Left panel shows a comparison between experimental and calculated values for the mole fraction standard Gibbs free energy of dimerization (∆ ) of the following transmembrane protein domains: ErbB1, EphA1, WALP23 and GpA. Simulation errors are estimated with

  tests involving complex systems are performed, including comparisons with experimental or atomistic simulation data. Here most of the system are considered validation. For a complete overview of the parametrization strategy used, see the Supplementary Notes, section B. The D B ΔGexp (kJ/mol) ΔGCG sim (kJ/mol) PMF (kJ/mol)

  , consistency was difficult to obtain given the limited bead types and sizes. To alleviate these problems, we undertook a rebalancing of all non-bonded interaction terms of the Martini model, including the addition of new beads and labels. The new version, called Martini 3, enables more accurate simulations of molecular systems in general. In this paper, we present the key features of Martini 3 combined with examples of new applications and improvements in relation to the previous Martini model.

  GROMACS 2016 and 2018)[START_REF] Abraham | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF][START_REF] Van Der Spoel | GROMACS: Fast, flexible, and free[END_REF] , VMD 1.9.4a12[START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF] , xmgrace (5.1.25) and MDAnalysis[START_REF] Gowers | MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations[END_REF] . The graphs were plotted using Excel 2016, xmgrace (5.1.25) and gnuplot (5.2). Figures were compiled using VMD 1.9.4a12 and Inkscape 1.1.

	(coupling time
	constant of 1.0 ps) and the Parrinello-Rahman barostat 86 (coupling time
	constant of 12.0 ps) were employed to maintain temperature and pressure,
	respectively. Except for equilibration runs, a time step of 20 fs was used for all
	systems. CG simulation settings are available as input files for GROMACS on
	the Martini portal http://cgmartini.nl. GROMACS 2016.x and 2018.x were
	used to run all the MD simulations 87,88 . For automated running and managing
	the Martini 3 simulations, we provide an updated version of Martinate 89,90 ,

described in Supplementary Codes, section H2. All the analysis were performed using gmx analysis tools (

ACKNOWLEDGMENTS

We thank all members of the S.J.M. group and also external users for testing Martini 3 in its open-beta version. In particular, we would like to thank C.F.E.

Schroer

, P.W.J.M. Frederix, W. Pezeshkian, M.N. Melo, H.I. Ingólfsson, M. Tsanai, M. König, P. A. Vainikka, T. Zijp, L. Gaifas, J.H. van der Woude, M. Espinoza Cangahuala, M. Scharte, J. Cruiming, L.M. van der Sleen, V. Verduijn, A.H. Beck Frederiksen, B. Schiøtt, M. Sikora, P. Schmalhorst, K. Pluhackova, C. Arnarez, C.A. López, E. Jefferys, and M.S.P. Sansom for their preliminary tests with a lot of different systems including aedamers, sugars, amino acids, deep eutectic solvents, lipids, peptides, and proteins. We also thank the Center for Information Technology of the University of Groningen for providing access to the Peregrine high-performance computing cluster. We acknowledge the National Computing Facilities Foundation (NCF) of The Netherlands Organization for Scientific Research (NWO), CSC -IT Center for Science Ltd (Espoo, Finland), and CINES (France) for providing computing time. Work in the S.J.M. group was supported by an ERC Advanced Grant "COMP-MICR-CROW-MEM". R.A. thanks The Netherlands Organisation for Scientific Research NWO (Graduate Programme Advanced Materials, No. 022.005.006) for financial support. L.M. acknowledges the Institut National de la Santé et de la Recherche Medicale (INSERM) and the Agence Nationale de la Recherche (ANR) for funding (grant no. ANR-17-CE11-0003) and GENCI-CINES for computing time (grant no. A0060710138). S.T. acknowledges the support from the European Commission via a Marie Skłodowska-Curie Actions individual fellowship (MicroMod-PSII, grant agreement 748895). M.J. thanks the Emil Aaltonen foundation for financial

support. I.V. thanks the Academy of Finland (Center of Excellence program (grant no. 307415)), Sigrid Juselius Foundation, the Helsinki Institute of Life Science fellow program, and the HFSP (research grant Ref.-No: RGP0059/2019). R.B. and J.D. were supported by the intramural research program of the NIDDK, NIH. Their work utilized the computational resources of the NIH HPC Biowulf cluster. (http://hpc.nih.gov). H.M.-S. acknowledges the Czech Science Foundation (19-19561S). J.B. acknowledges funding from the TOP grant of S.J.M. (NWO) and the EPSRC program grant EP/P021123/1. Work in D.P.T.'s group is supported by the Natural Sciences and Engineering Research Council (Canada) and Compute Canada, funded by the Canada Foundation for Innovation. D.P.T. acknowledges further support from the Canada Research Chairs program.

DATA AVAILABILITY

Force-field parameters and procedures (e.g. tutorials) are publicly available at http://cgmartini.nl. Simulation Data (e.g. trajectories) supporting the results of this paper are available from the corresponding authors upon reasonable request.

CODE AVAILABILITY

Martinize2 (for which the manuscript is in preparation) and Martinate codes used in this work are publicly available at https://github.com/marrink-lab/. For a more detailed information, see Supplementary Codes, section H.

AUTHOR CONTRIBUTIONS

P.C.T.S. and S.J.M. conceived the project with suggestions from R.A., A.H.V., J.B., and S.T.; P.C.T.S. generated and optimized all bead parameters; P.C.T.S., R.A., and J.B. generated the topology and bonded parameters of all CG models with suggestions from S.T. and I.F.; P.C.T.S., R.A., A.H.V., and F.G. performed the simulations and analysis involving transfer free energies, solvent and polymer properties; P.C.T.S., S.T., J.B., and J.M. performed the simulations and analysis involving lipid bilayers; P.C.T.S., I.F., and R.A. performed the simulations and analysis involving nucleobases; P.C.T.S., I.P., and A.H.V. generated the models, performed the simulations and analysis involving aedamers; P.C.T.S., and F.G. generated the models, performed the simulations and analysis involving ionic liquids and ionic water solutions; R.A. generated the models, performed the simulations and analysis involving bulk heterojunctions, with suggestions from L.M. regarding the fullerene model; P.C.T.S., J.B., H.A., R.A., B.M.H.B., S.T., J.M., V.N., X.P., M.J., H.M.K., J.D., V.C., and H.M.-S. performed the simulations and analysis involving amino acids, peptides and proteins; J.B., T.W., P.K., and S.T. developed some tools and scripts used to generate the CG models and to run the MD simulations;

COMPETING INTERESTS

The authors declare no competing interests.

ONLINE METHODS

1) CG models

CG mappings of small molecules were initially inspired by the standard Martini 2 models, when they were available. Due to the well-balanced properties of the regular (R), small (S), and tiny (T) beads in Martini 3, the CG models now follow more specific rules for mapping. For instance, over-representing 3-to-1 or 2-to-1 fragments by the usage of R-beads is always avoided. Aromatic rings without substituents are composed of T-beads and, in case of substituents, S-beads are used. Aliphatic rings without substituents are usually based on S-beads, which better reproduce their molecular shape.

More technical details about the mapping rules and bead types used are given in the Supplementary Notes, sections C1 and C3. As in the previous version of Martini [START_REF] Uusitalo | Coarse-Grained Force Field: Extension to DNA[END_REF][START_REF] Uusitalo | Martini Coarse-Grained Force Field: Extension to RNA[END_REF][START_REF] López | Martini Coarse-Grained Force Field: Extension to Carbohydrates[END_REF][START_REF] Monticelli | The MARTINI Coarse-Grained Force Field: Extension to Proteins[END_REF][START_REF] Grunewald | Transferable MARTINI Model of Poly(ethylene Oxide)[END_REF] , bonded parameters are based on atomistic simulations or high-resolution experimental data. The main difference in Martini 3 lies in the protocol to derive bond lengths, which are now based on matching overall volume and shape of the molecules (Supplementary Notes, section C2). In this spirit, the bonded parameters of the protein models were also slightly modified from the standard Martini 2.2 values [START_REF] Monticelli | The MARTINI Coarse-Grained Force Field: Extension to Proteins[END_REF][START_REF] De Jong | Improved Parameters for the Martini Coarse-Grained Protein Force Field[END_REF] , including the addition of side chain corrections [START_REF] Herzog | Improved Side Chain Dynamics in MARTINI Simulations of Protein-Lipid Interfaces[END_REF] , based on experimental reference structures. Backbone bead types do not depend on the secondary structure anymore, but are now represented by P2 beads, except for proline (SP1a), alanine (SP2, with an additional bead for the side chain) and glycine (SP1).

Adapted versions of Gō-like models [START_REF] Poma | Combining the MARTINI and Structure-Based Coarse-Grained Approaches for the Molecular Dynamics Studies of Conformational Transitions in Proteins[END_REF] or Elastic Networks [START_REF] Periole | Combining an Elastic Network With a Coarse-Grained Molecular Force Field: Structure, Dynamics, and Intermolecular Recognition[END_REF] were used to maintain the tertiary protein structure. All CG protein models were constructed using Martinize2, described in Supplementary Codes, section H1. Lipid Supplementary Notes, section D provide details of the specific systems and methods used in the tests performed to parametrize and validate the new Martini 3 LJ parameters. Further information on research design is available in the Life Sciences Reporting Summary linked to this article.