A Cecile

J W Janssens
email: a.janssens@erasmusmc.nl

John P A Ioannidis

Cornelia M Van Duijn

Julian Little

Muin J Khoury

Strengthening the Reporting of Genetic Risk Prediction Studies: The GRIPS Statement

Keywords:

Introduction

The recent successes of genome-wide association studies and the promises of whole genome sequencing fuel interest in the translation of this new wave of basic genetic knowledge to health care practice. Knowledge about genetic risk factors may be used to target diagnostic, preventive, and therapeutic interventions for complex disorders based on a person's genetic risk, or to complement existing risk models based on classical nongenetic factors such as the Framingham risk score for cardiovascular disease. Implementation of genetic risk prediction in health care requires a series of studies that encompass all phases of translational research [START_REF] Khoury | The continuum of translation research in genomic medicine: how can we accelerate the appropriate integration of human genome discoveries into health care and disease prevention?[END_REF][START_REF] Hlatky | Criteria for evaluation of novel markers of cardiovascular risk: a scientific statement from the American Heart Association[END_REF], starting with a comprehensive evaluation of genetic risk prediction.

With increasing numbers of discovered genetic markers that can be used in future genetic risk prediction studies, it is crucial to enhance the quality of the reporting of these studies, since valid interpretation could be compromised by the lack of reporting of key information. Information that is often missing includes details in the description of how the study was designed and conducted (e.g., how genetic variants were selected and coded, how risk models or genetic risk scores were constructed, and how risk categories were chosen), or how the results should be interpreted. An appropriate assessment of the study's strengths and weaknesses is not possible without this information. There is ample evidence that prediction research often suffers from poor design and bias, and these may also have an impact on the results of the studies and on models of disease outcomes based on these studies [START_REF] Kyzas | Quality of reporting of cancer prognostic marker studies: association with reported prognostic effect[END_REF][START_REF] Kyzas | Selective reporting biases in cancer prognostic factor studies[END_REF][START_REF] Mcshane | REporting recommendations for tumor MARKer prognostic studies (REMARK)[END_REF].

Although most prognostic studies published to date claim significant results [START_REF] Kyzas | Almost all articles on cancer prognostic markers report statistically significant results[END_REF][START_REF] Tzoulaki | Assessment of claims of improved prediction beyond the Framingham risk score[END_REF], very few translate to clinically useful applications. Just as for observational epidemiological studies [START_REF] Von Elm | The scandal of poor epidemiological research[END_REF], poor reporting complicates the use of the specific study for research, clinical, or public health purposes and hampers the synthesis of evidence across studies.

Reporting guidelines have been published for various research designs [START_REF] Simera | A catalogue of reporting guidelines for health research[END_REF], and these contain many items that are also relevant to genetic risk prediction studies. In particular, the guidelines for genetic association studies (STREGA) have relevant items on the assessment of genetic variants, and the guidelines for observational studies (STROBE) have relevant items about the reporting of study design. The guidelines for diagnostic studies (STARD) and those for tumor marker prognostic studies (REMARK) include relevant items about test evaluation; the REMARK guidelines also have relevant items about risk prediction [START_REF] Von Elm | The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies[END_REF][START_REF] Little | STrengthening the REporting of Genetic Association Studies (STREGA): an extension of the STROBE statement[END_REF][START_REF] Bossuyt | Towards complete and accurate reporting of studies of diagnostic accuracy: the STARD initiative[END_REF][START_REF] Mcshane | Reporting recommendations for tumor marker prognostic studies[END_REF]. However, none of these guidelines are fully suited to genetic risk prediction studies, an emerging field of investigation with specific methodological issues that need to be addressed, such as the handling of large numbers of genetic variants (from 10s to 10,000s) and flexibility in handling such large numbers in analyses. We organized a two-day workshop with an international group of risk prediction researchers, epidemiologists, geneticists, methodologists, statisticians, and journal editors to develop recommendations for the reporting of Genetic RIsk Prediction Studies (GRIPS).

Genetic Risk Prediction Studies

Genetic risk prediction studies typically develop or validate models that predict the risk of disease, but they are also being investigated for use in predicting prognostic outcome, treatment response, or treatment-related harms. Risk prediction models are statistical algorithms, which may be simple genetic risk scores (e.g., risk allele counts), may be based on regression analyses (e.g., weighted risk scores or predicted risks), or may be based on more complex analytic approaches such as support vector machine learning or classification trees. The risk models may be based on genetic variants only, or include both genetic and nongenetic risk factors [START_REF] Janssens | Genome-based prediction of common diseases: methodological considerations for future research[END_REF].

Aims and Use of the GRIPS Statement

The 25 items of the GRIPS statement are intended to maximize the transparency, quality, and completeness of reporting on research methodology and findings in a particular study. It is important to emphasize that these recommendations are guidelines only for how to report research and do not prescribe how to perform genetic risk prediction studies. The guidelines do not support or oppose the choice of any particular study design or method, e.g., the guidelines recommend that the study population should be described, but do not specify which population is preferred in a particular study.

The intended audience for the reporting guidelines is broad and includes epidemiologists, geneticists, statisticians, clinician scientists, and laboratory-based investigators who undertake genetic risk prediction studies, as well as journal editors and reviewers who have to appraise the design, conduct and analysis of such studies.

In addition, it includes -users‖ of such studies who wish to understand the basic premise, design, and limitations of genetic prediction studies in order to interpret the results for their potential application in health care. These guidelines are also intended to ensure that essential data from future genetic risk prediction studies are presented in standardized form, which will facilitate information synthesis as part of systematic reviews and meta-analyses.

Items presented in the checklist are relevant for a wide array of risk prediction studies, because GRIPS focuses on the main aspects of the design and analysis of risk prediction studies. GRIPS does not address randomized trials that may be performed to test risk models, nor does it specifically address decision analyses, costeffectiveness analyses, assessment of health care needs, or assessment of barriers to health care implementation [START_REF] Khoury | The emergence of translational epidemiology: from scientific discovery to population health impact[END_REF]. Once the performance of a risk model has been established, these next steps toward implementation require further evaluation [START_REF] Von Elm | The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies[END_REF][START_REF] Moons | Prognosis and prognostic research: application and impact of prognostic models in clinical practice[END_REF]. For the reporting of these studies, which go beyond the assessment of genetic risk models as such, additional requirements apply. However, proper documentation of genetic predictive research according to GRIPS might facilitate the translation of research findings into clinical and public health practice.

Development of the GRIPS Statement

The GRIPS statement was developed by a multidisciplinary panel of 25 risk prediction researchers, epidemiologists, geneticists, methodologists, statisticians, and journal editors, seven of whom were also part of the STREGA initiative [START_REF] Little | STrengthening the REporting of Genetic Association Studies (STREGA): an extension of the STROBE statement[END_REF]. They attended a two-day meeting in Atlanta, Georgia (US) in December 2009 that was sponsored by the US Centers for Disease Control and Prevention on behalf of the Human Genome Epidemiology Network (HuGENet) [START_REF] Khoury | The Human Genome Epidemiology Network[END_REF]. Participants discussed a draft version of the guidelines that was prepared and distributed before the meeting.

This draft version was developed on the basis of existing reporting guidelines, namely STREGA [START_REF] Little | STrengthening the REporting of Genetic Association Studies (STREGA): an extension of the STROBE statement[END_REF], REMARK [START_REF] Mcshane | Reporting recommendations for tumor marker prognostic studies[END_REF], and STARD [START_REF] Bossuyt | Towards complete and accurate reporting of studies of diagnostic accuracy: the STARD initiative[END_REF]. These were selected out of all available guidelines (see http://www.equator-network.org) because of their focus on observational study designs and genetic factors (STREGA), prediction models (REMARK), and test evaluation (REMARK and STARD). During the meeting, methodological issues pertinent to risk prediction studies were addressed in presentations. Workshop participants were asked to change, combine, or delete proposed items and add additional items if necessary. Participants had extensive post-meeting electronic correspondence. To harmonize our recommendations for genetic risk prediction studies with previous guidelines, we chose the same wording for the items wherever possible. Finally, we tried to create consistency with previous guidelines for the evaluation of risk prediction studies of cardiovascular diseases and cancer [START_REF] Hlatky | Criteria for evaluation of novel markers of cardiovascular risk: a scientific statement from the American Heart Association[END_REF][START_REF] Freedman | Cancer risk prediction models: a workshop on development, evaluation, and application[END_REF]. The final version of the checklist is presented in Table 1.

The GRIPS Explanation and Elaboration Article

Accompanying this GRIPS statement, an Explanation and Elaboration document has been written (see Text S1), modeled after those developed for other reporting guidelines [START_REF] Altman | The revised CONSORT statement for reporting randomized trials: explanation and elaboration[END_REF][START_REF] Bossuyt | The STARD statement for reporting studies of diagnostic accuracy: explanation and elaboration[END_REF][START_REF] Liberati | The PRISMA statement for reporting systematic reviews and meta-analyses of studies that evaluate health care interventions: explanation and elaboration[END_REF][START_REF] Vandenbroucke | Strengthening the Reporting of Observational Studies in Epidemiology (STROBE): explanation and elaboration[END_REF]. The Explanation and Elaboration document illustrates each item with at least one published example that we consider transparent in reporting, explains the rationale for its inclusion in the checklist, and presents details of the items that need to be addressed to ensure transparent reporting. The Explanation and Elaboration document was produced after the meeting. The document was prepared by a small subgroup and shared with all workshop participants for additional revisions and final approval.

Concluding Remarks and Future Directions

High-quality reporting reveals the strengths and weaknesses of empirical studies, facilitates the interpretation of the scientific and health care relevance of the resultsespecially within the framework of systematic reviews and meta-analyses-and helps build a solid evidence base for moving genomic discoveries into applications in health care practice. The GRIPS guidelines were developed to improve the transparency, quality and completeness of the reporting of genetic risk prediction studies. As outlined in the introduction, GRIPS does not prescribe how studies should be designed, conducted, or analyzed, and therefore the guidelines should not be used to assess the quality of empirical studies [START_REF] Vandenbroucke | and CONSORT: for whom does the guideline toll?[END_REF]. The guidelines should be used only to check whether all essential items are adequately reported.

Finally, the methodology for designing and assessing genetic risk prediction models is still developing. For example, newer measures of reclassification were first introduced in 2007 [START_REF] Cook | Use and misuse of the receiver operating characteristic curve in risk prediction[END_REF], and several alternative reclassification measures have been proposed [START_REF] Pencina | Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond[END_REF]. Which measures to apply and when to use measures of reclassification are still subject to ongoing evaluation and discussion [START_REF] Janssens | Assessment of improved prediction beyond traditional risk factors: when does a difference make a difference?[END_REF]. Furthermore, alternative strategies for constructing risk models other than simple regression analyses are being explored, and these may add increased complexity to the reporting. In formulating the items of the GRIPS statement, these methodological advances were anticipated. It is for this reason that the GRIPS statement recommends how a study should be reported and not how a study should be conducted or analyzed. Therefore, methodological and analytical developments will not immediately impact the validity and relevance of the items, but the GRIPS statement will be updated when this is warranted by essential new developments in the construction and evaluation of genetic risk models. * Marked items should be reported for every population in the study.

Table 1 .

 1 Reporting recommendations for evaluations of risk prediction models that include genetic variants. Report distributions of predicted risks and/or risk scores. Assessment 18 Report measures of model fit and predictive ability, and any other performance measures, if pertinent. Validation 19 Report any validation of the risk model(s). Other analyses 20 Present results of any subgroup, interaction, or exploratory analyses, whenever pertinent. DISCUSSION Limitations 21 Discuss limitations and assumptions of the study, particularly those concerning study design, selection of participants, and measurements and analyses, and discuss their impact on the results of the study. Interpretation 22 Give an overall interpretation of results considering objectives, limitations, multiplicity of analyses, results from similar studies, and other relevant evidence. Generalizability 23 Discuss the generalizability and, if pertinent, the health care relevance of the study results. OTHER Supplementary information 24 State whether databases for the analyzed data, risk models, and/or protocols are or will become publicly available and if so, how they can be accessed. Funding 25 Give the source of funding and the role of the funders for the present study. State whether there are any conflicts of interest.

	TITLE & ABSTRACT
		1	(a) Identify the article as a study of risk prediction using genetic factors. (b) Use
			recommended keywords in the abstract: genetic or genomic, risk, prediction.
	INTRODUCTION
	Background and	2	Explain the scientific background and rationale for the prediction study.
	rationale
	Objectives	3	Specify the study objectives and state the specific model(s) that is/are investigated. State
			if the study concerns the development of the model(s), a validation effort, or both.
	METHODS
	Study design and	4* Specify the key elements of the study design and describe the setting, locations, and
	setting		relevant dates, including periods of recruitment, follow-up, and data collection.
	Participants	5* Describe eligibility criteria for participants, and sources and methods of selection of

Acknowledgments

The members of the GRIPS group are: A. Cecile J. W. Janssens, John P. A. Ioannidis, Sara Bedrosian, Paolo Boffetta, Siobhan M. Dolan, Nicole Dowling, Isabel Fortier, Andrew N. Freedman, Jeremy M. Grimshaw, Jeffrey Gulcher, Marta Gwinn, Mark A. Hlatky, Holly Janes, Peter Kraft, Stephanie Melillo, Christopher J. O'Donnell, Michael J. Pencina, David Ransohoff, Sheri D. Schully, Daniela Seminara, Deborah M. Winn, Caroline F. Wright, Cornelia M. van Duijn, Julian Little, and Muin J. Khoury.

Workshop was sponsored by the Centers for Disease Control and Prevention on behalf of the Human Genome Epidemiology Network (HuGENet). The findings and conclusions in this report are those of the authors and do not necessarily reflect the views of the Department of Health and Human Services. A. Cecile J.W. Janssens is financially supported by grants from the Erasmus University Medical Center Rotterdam, the Center for Medical Systems Biology in the framework of the Netherlands Genomics Initiative (NGI) and the VIDI grant of the Netherlands Organisation for Scientific Research (NWO). John P.A. Ioannidis: Tufts CTSI is supported by the National Institutes of Health/National Center for Research Resources (UL1 RR025752). Opinions in this paper are those of the authors and do not necessarily represent the official position or policies of the Tufts CTSI. Julian Little holds a Canada Research Chair in Human Genome Epidemiology. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Summary Points

The rapid and continuing progress in gene discovery for complex diseases is fueling interest in the potential application of genetic risk models for clinical and public health practice.

The number of studies assessing the predictive ability is steadily increasing, but the quality and completeness of reporting varies.

A multidisciplinary workshop sponsored by the Human Genome Epidemiology Network developed a checklist of 25 items recommended for strengthening the reporting of Genetic RIsk Prediction Studies (GRIPS), building on the principles established by prior reporting guidelines.

These recommendations aim to enhance the transparency of study reporting, and thereby to improve the synthesis and application of information from multiple studies that might differ in design, conduct, or analysis.

A detailed Explanation and Elaboration document is published as supporting information (Text S1).

Supporting Information

Text S1

Strengthening the Reporting of Genetic RIsk Prediction Studies (GRIPS): Explanation and Elaboration.


Yong Sheng

Dongmin Yang

Yuanqiang Tan

Jianqiao Ye
email: j.ye@leeds.ac.uk

Microstructure effects on transverse cracking in composite laminae by DEM

Keywords: Discrete element method, Transverse cracking, Debonding, Matrix cracking, Microstructure, Composite laminae

A particle discrete element method (DEM) was employed to simulate transverse cracking in laminated fiber reinforced composites. The microstructure of the laminates was modeled by a DEM model using different mechanical constitutive laws and materials parameters for different constituents, i.e., fiber, matrix and fiber/matrix interface. Rectangular, hexagonal and random fiber distributions were simulated to study the effect of fiber distribution on the transverse cracking. The initiation and dynamic propagation of transverse cracking and interfacial debonding were all captured by the DEM simulation, which showed similar patterns to those observed from experiments. The effect of fiber volume fraction was also studied for laminae with randomly distributed fibers. It was found that the distribution and volume fraction of fibers affected not only the transverse cracking path, but also the behavior of matrix plastic deformation and fiber/matrix interface yielding in the material.

Introduction

The inherent multi-phased microstructure of fiber reinforced composites, e.g. fiber volume fraction, fiber direction and interfacial strength, results in complex failure modes dominated by fracture of constituent elements. The damage evolution and failure mechanism of composites are obviously far more complex than those of monolithic materials due to the interactions of constituents. The complicated nature of failure/damage modes and physical composition of composite materials make the prediction and evaluation of progressive damage a challenging task. It is very important for material design and structural optimization in composite manufacturing to have the predictive capability on the dynamic failure/damage process. For instance, the properties of a laminate depend very much on the properties of its lamina and ply interface; the behavior of each lamina is, in turn, governed by its constituents, i.e. the properties of the fibers, the surrounding matrix, the interface and the relative amount of fibers and matrix in the lamina. When a composite is under external loading, damage, such as delamination, debonding, matrix cracking and fiber breakage, may occur. The accumulation of the damage during the loading process may eventually cause collapse of the whole material. Therefore a study of the damage or failure modes at micro scale with the consideration of microstructure is a more direct approach to predict the mechanical behavior of the composites.

Traditional micromechanical models focused mainly on the analysis of stress field [START_REF] Asp | Prediction of matrix-initiated transverse failure in polymer composites[END_REF][START_REF] Asp | Effects of fiber and interphase on matrix-initiated transverse failure in polymer composites[END_REF][START_REF] Trias | Random models versus periodic models for fibre reinforced composites[END_REF], or predicting crack propagation using a prescribed failure criterion for composites with small numbers of fibers [START_REF] Zhang | Effects of fibre debonding and sliding on the fracture behaviour of fibre-reinforced composites[END_REF][START_REF] Rao | Micro-mechanical modeling of machining of FRP composites -Cutting force analysis[END_REF][START_REF] Vejen | Transverse crack growth in glass/epoxy composites with exactly positioned long fibres[END_REF]. The ability to cope with inhomogeneity or/and disorder is essential in order to predict the damage process of a composite from micro cracking to catastrophic failure. Transverse cracking, which normally is formed by matrix cracking and fiber/matrix interfacial debonding, is a familiar phenomenon in a lamina or laminated composite when it is subjected to transverse tension. Finite element method (FEM) has been widely used in this analysis on the basis of imposing various failure criteria [START_REF] Berthelot | Transverse cracking of cross-ply laminates: Part 2. Progressive widthwise cracking[END_REF][START_REF] Sjögren | The effects of matrix and interface on damage in GRP cross-ply laminates[END_REF][START_REF] Okabe | Numerical modeling of progressive damage in fiber reinforced plastic cross-ply laminates[END_REF]. However, many of these criteria do not have a rigorous physical basis that can be related to the microstructure of composite materials. Furthermore, the stress field that is calculated to meet the failure criteria of a composite is normally obtained by assuming that the material is not damaged. Thus stress redistribution in the material arising from, e.g., strain softening induced by localized damage cannot be taken into account. Some of the FEM models may agree well with experimental results, but it is very difficult for them to explain how the damage starts and develops in a composite during the entire loading process. Therefore, it is very challenging for a conventional continuum mechanics based numerical model to predict the whole dynamical damage/failure process, from matrix cracking initiation to final laminae/laminates collapse.

A fiber reinforced composites can be imaged as discrete fibers being inserted into a continuous matrix and bonded together by the interfaces between them. The discrete fibers interact with each other through the matrix between them. The analysis of the composite can be carried out through the analysis of the interaction of the discrete fibers. Since a large number of fibers are used in a fiber reinforced composite material, the interaction of these fibers was conventionally studied by the Monte-Carlo method [START_REF] Okabe | Characterization of tensile damage and strength in GFRP cross-ply laminates[END_REF][START_REF] Mishnaevsky | Statistical modelling of compression and fatigue damage of unidirectional fiber reinforced composites[END_REF]. However, the Monte Carlo method cannot take into account the matrix effect, and the numerical results can hardly interpret the experimental investigations. On the other hand, discrete element method (DEM) [START_REF] Cundall | A discrete mumerical model for granular assemblies[END_REF] combined the statistical analysis with the concept of physical bonds from Molecular Dynamics (MD), and offered an intrinsic advantage over other numerical methods in modeling microscopic and multiphased materials. DEM was proposed by Cundall in the context of rock mechanics [START_REF] Cundall | A discrete mumerical model for granular assemblies[END_REF] and has been implemented in many other fields, such as geomaterials [START_REF] Hunt | Modelling the Kaiser effect and deformation rate analysis in sandstone using the discrete element method[END_REF], granular materials [START_REF] Sheng | Numerical studies of uniaxial powder compaction process by 3D DEM[END_REF], concrete [START_REF] Frédéric | Discrete element modeling of concrete submitted to dynamic loading at high strain rates[END_REF], ceramics [START_REF] Tan | Study of polycrystalline Al2O3 machining cracks using discrete element method[END_REF][START_REF] Tan | Discrete element method (DEM) modeling of fracture and damage in the machining process of polycrystalline SiC[END_REF], particle-fluid [START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows: Computational issues[END_REF] and particle-gas two phase flow [START_REF] Chu | CFD-DEM modelling of multiphase flow in dense medium cyclones[END_REF]. Kim et al. [20,[START_REF] Kim | Numerical fracture analysis on the specimen size dependency of asphalt concrete using a cohesive softening model[END_REF] used DEM to simulate the damage of asphalt concrete which is a typical example of two-phase composite materials. Yang et al [START_REF] Yang | Discrete element modeling of the microbond test of fiber reinforced composite[END_REF] used DEM to simulate the dynamic process of microbond test of single fiber reinforced composite with the consideration of matrix cracking as well as fiber/matrix interfacial debonding. Donze et al. [START_REF] Magnier | Numerical simulations of impacts using a discrete element method[END_REF][START_REF] Shiu | Penetration prediction of missiles with different nose shapes by the discrete element numerical approach[END_REF] simulated the impacts of rigid spherical nose shaped missiles on concrete beams with reinforcements. The existing studies concluded that DEM could provide unique insight into damage evolution at both micro and macro scales in the process of the dynamic characterization of fracturing. Van Mier et al. [START_REF] Van Mier | Influence of microstructure of concrete on size/scale effects in tensile fracture[END_REF][START_REF] Van Mier | Some Notes on Microcracking, Softening, Localization, and Size Effects[END_REF] used a lattice beam model, which is equivalent to the DEM model using regular particle arrangements, to study the fracture behavior of concrete from micro-cracking initiation and macro-crack failure in both two and three dimensions. Wittel et al. [START_REF] Wittel | Discrete element simulation of transverse cracking during the pyrolysis of carbon fibre reinforced plastics to carbon/carbon composites[END_REF][START_REF] Wittel | A study of transverse ply cracking using a discrete element method[END_REF] has applied a lattice spring model, which is a simple form of DEM, to simulate ply transverse cracking of cross-ply laminates, in which the nodes represent fibers and the springs with random breaking thresholds represent the disordered matrix. It has been confirmed by the existing studies that DEM is particularly suitable to simulate dynamic instability in crack propagation and the collective behaviour of many interacting cracks. However, the spring models used for the matrix and interface by Wittel et al. [START_REF] Wittel | Discrete element simulation of transverse cracking during the pyrolysis of carbon fibre reinforced plastics to carbon/carbon composites[END_REF][START_REF] Wittel | A study of transverse ply cracking using a discrete element method[END_REF] are too simple to represent the mechanical properties of real materials and interfaces. Matrix cracking cannot be correctly simulated by this model and, only regular distribution of fibers with a fixed volume fraction can be considered.

In the research presented in this paper, DEM was used to simulate transverse cracking in laminae, including onset and propagation of the matrix cracking and fiber/matrix interface debonding. Composite laminae with regular distributions of fibers were modeled first and the numerical results were compared with the existing experimental results to validate the developed DEM model. The model was then modified to include random distribution of fibers with different fiber volume fractions. Based on the numerical results from the DEM modelling, the effect of microstructure of the composite laminae on the onset and the propagation path of transverse cracking, and the resulting residual damage within the material can be assessed. The DEM model proposed in this paper can be used for the analysis of more complex composite systems, such as cross-ply laminates.

Discrete Element Method (DEM)

The particle discrete element method used in this research assumes that the particle elements are rigid spheroids (3D) or discs (2D), and can overlap or detach from each other. The contact forces between any two particles are determined from the overlap and relative movements of the particle pair according to a specified force-displacement law. In 2D DEM, the motion of a particle over a time step ∆t is governed by the Newton's second law as below [START_REF] Cundall | A discrete mumerical model for granular assemblies[END_REF][START_REF] Sheng | Numerical studies of uniaxial powder compaction process by 3D DEM[END_REF]:

Translational motion (1)
Rotational motion [START_REF] Asp | Effects of fiber and interphase on matrix-initiated transverse failure in polymer composites[END_REF] where i (= 1, 2) denotes, respectively, the x-and y-co-ordinate directions; F i is the resultant force of the particle; v i is the translational velocity; m is the mass of the particle;. g i is the body force acceleration vector (e.g., gravity loading); M 3 is the out of balance moment referred to the out-of-plane axis; w 3 is the rotational velocity about the out-of-plane axis; I is the rotational inertia of the particle; t is time and ∆ is the increment. Damping, e.g., local damping or viscous damping, can be added in the DEM model to dissipate the kinetic energy together with particles' frictional sliding, such that a steady-state solution can be obtained more efficiently [START_REF]PFC2D (particle flow code in 2-dimensions), Version 3.10[END_REF]. Equations ( 1) and ( 2) are solved by a finite difference scheme. Both the specified force-displacement law and Newton's second motion law are used in the calculation cycle of the discrete element method.

DEM allows particles to be bonded together at contacts and to separate when the bond strength or energy is exceeded. Therefore it can simulate the motion of individual particles and also the behavior of a bulk material which is formed by assembling many particles through bonds at contacts with different constitutive laws. In a DEM model of bulk material, elementary micro scale particles are assembled to form the bulk material with macroscopic continuum behavior determined only by the dynamic interaction of particles [START_REF]PFC2D (particle flow code in 2-dimensions), Version 3.10[END_REF][START_REF] Potyondy | A bonded-particle model for rock[END_REF]. Unlike the traditional solution using the strain and stress relations, contact properties are the dominating parameters in a DEM solution, combined with size and shape of the particles. Subject to external loading, when the strength or the fracture energy of a bond between particles is exceeded, flow and disaggregation of the particle assembly occur and the bond starts breaking. Consequently, cracks form naturally at micro scale. Hence, damage modes and their interaction emanate as the process of debonding of particles. The way that DEM discretizes a material domain gives the most significant advantage over the traditional continuum methodologies, as dynamic material behavior of composites, crack tip singularities and crack formulation criteria can all be avoided due to the naturally discontinuous representation of composite materials via particle assemblies.

In a DEM model for composite materials, damage initiation and propagation will take place at the particle contacts under the applied load. There is no requirement for additional failure criteria and the damage of different constituents can be taken into account simultaneously. The damage interaction will also be addressed through the dynamic motion of particles and contacts. The simulation results of a DEM model are naturally discrete and collected at microscopic level. The statistical analysis of the DEM results can also reveal material and structural behavior at macro scale by studying the bonds that have been broken or has yielded.

DEM Simulation of Laminae under Transverse Tension

Constitutive models

A DEM model of a solid material is usually constructed by packing discrete particles into an assembly through rectangular, hexagonal or random arrangement. Random arrangement is the most general approach in which the disperse particles are compacted or expanded to reach a dense contacting network, as shown in Fig. 1.

Since a lamina is an inhomogeneous material which includes three phases: fiber, matrix and fiber/matrix interface, individual DEM models for each of the three phases are developed and virtual material tests (e.g., tensile and compressive tests) are carried out on all the phases, respectively, before the construction of the whole lamina model. The fiber material is represented by a parallel bond model whose constitutive behavior is illustrated in Fig. 2. A parallel bond can be regarded as a set of elastic springs with constant normal and shear stiffness, uniformly distributed over either a circular or rectangular cross-section lying on the contact plane and centered at the contact point. The bond can transmit both forces and moments. A contact softening model which is similar to the cohesive zone model (CZM) [START_REF] Hu | Stable numerical simulations of propagations of complex damages in composite structures under transverse loads[END_REF] in continuum mechanics is used to account for the matrix and fiber/matrix interface. The only difference between these two models is the unloading and reloading paths after yielding, as shown in Fig. 3. The contact softening model describes plastic deformation by linearly softening the bond after the contact force reaches the bond strength. In both tension and shearing, the bond strength decreases to zero when the plastic displacement reaches the maximum plastic displacement U pmax that is related to the fracture energy release rate G. The interfacial crack may be classified as mode I, mode II or mixed mode according to the stress field at the crack tip. In order to automatically capture the three fracture models, the maximum plastic displacement U pmax is kept constant, while the normal and shear strengths of the bond are defined individually. Hence, in a two dimensional system, the fracture energy release rate for mode I and mode II are, respectively:

(

In 2D DEM, the contact stresses at the bond are taken as the average stress between the two contacting particles [START_REF]PFC2D (particle flow code in 2-dimensions), Version 3.10[END_REF], as calculated below:

(

Where, R 1 and R 2 are the radii of the two contacting particles, respectively; δ is the dimension in the out-of-plane direction. k is in place of n or s that indicates normal or shear direction, respectively. Under tension or shearing , the bond starts to yield when the contacting force exceeds the strength and eventually breaks when the plastic displacement reaches the maximum plastic displacement U pmax . In the situation of mixed mode fracture, the bond strength F max (in force unit), is calculated from the two strength parameters (i.e. F c n and F c s , also in force unit) as a function of the current direction of the contact force. It is assumed that the contact strength varies linearly with the angle, α:

(7)
Where α is the angle between the directions of the contact force and the line segment connecting the centers of the two particles. The yielding of the bond in tension is determined by comparing the resultant contact force, i.e.

(

) 8
with the contact strength. The bond yields if the contact force is larger than the contact strength:

The maximum plastic displacement U pmax is calculated from the displacements in the normal and shear directions. In each computation step, the increments of normal and shear plastic deformation, which are denoted by ∆, are subjected to a consistence condition of:

When the plastic displacement reaches U pmax ,

The fracture energy release rate can be calculated as: [START_REF] Hunt | Modelling the Kaiser effect and deformation rate analysis in sandstone using the discrete element method[END_REF] In which σ n and σ s are the normal and shear stresses of the bond when yield occurs, i.e., equations ( 8) and ( 9) are satisfied. For the mixed mode, the fracture energy release rate is somewhere between the rates of the two single fracture modes.

Convergence study

Since displacement loading is used in this research, including the fiber, matrix and the laminae tensile tests shown in the following sections, the influence of loading velocities on the solutions is studied by comparing the stress-strain curves obtained from a parallel bond model subject to different tensile loading velocities, as shown in Fig. 4. It is found that for all loading velocities below 0.1m/s, the tests lead to almost identical results. Therefore, v=0.1m/s is chosen as the loading velocity in all subsequent cases to achieve the quasi-static conditions. After the convergence study on velocity, a Double Cantilever Beam (DCB) virtual test [START_REF] Meo | Delamination modelling in a double cantilever beam[END_REF][START_REF] Li | Use of mode-I cohesive-zone models to describe the fracture of an adhesively-bonded polymer-matrix composite[END_REF], which has been widely applied to investigate interface strength, is conducted by the 2D DEM model with square arrangement of particles. The test is to assess the effect of particle size on the failure of the contact softening model. The model consists of two identical material layers [START_REF] Ouyang | Local Damage Evolution of Double Cantilever Beam Specimens During Crack Initiation Process: A Natural Boundary Condition Based Method[END_REF], as shown in Fig. 5, each of which is further divided into 2 plies. The contacts between the particles within each layer are modeled by the parallel bond, while the interface between the two layers is represented by the contact softening model. The length of the plies is 45mm and the thickness is 3mm. An initial interfacial crack of 20mm length from the left end of the beam is generated by removing the interfacial bonds between the two layers. The right end of the beam is fixed. The loading is applied by assigning a constant separation velocity in the vertical direction to the particles at the left ends of the two layers, which results in a variable opening end displacement. Different particle sizes are used in this simulation to examine their effect on the interfacial fracture and the results are illustrated in Fig. 6.

U n U n p U n pmax F n F c n U s U s p U s pmax F c s F s (
It is found that the model shows good convergence and numerical results agree very well with the theoretical solution by the J-integral theory in [START_REF] Ouyang | Local Damage Evolution of Double Cantilever Beam Specimens During Crack Initiation Process: A Natural Boundary Condition Based Method[END_REF], though the results for R=0.75mm are slightly distanced from the group, which is due to the fact that for this case only one ply of particles is used in each layer. The good convergence of the contact softening model indicates that the influence of particle size on the model is not significant as long as the particle size is sufficiently small and the model can be adopted to simulate the fracture of matrix and fiber/matrix interface.

DEM simulations

Glass fiber reinforced composites are taken as examples in this paper, where the glass fiber and the epoxy matrix are both considered as linear isotropic materials. The matrix, the fiber and then the composites are all generated using the particle assembly model [START_REF] Tan | Study of polycrystalline Al2O3 machining cracks using discrete element method[END_REF][START_REF]PFC2D (particle flow code in 2-dimensions), Version 3.10[END_REF]. Particles with a radius between R min and R max are randomly generated according to the normal distribution law in a cell. Since the particles are randomly packed, it is not realistic to derive directly a rigorous and theoretical relationship between the micro properties and the macro mechanical properties, such as between particle stiffness, contact stiffness and contact strength and the elastic modulus, Poisson's ratio and tensile strength. However, each of the macro mechanical properties is related to one or several of the micro parameters. For instance, the elastic modulus is only related to the particle and contact stiffness, but not the contact strength and particle size. Based on this observation, virtual simulation tests are usually carried out to determine the micro properties in order to obtain the required macro mechanical properties. Virtual tensile tests on the DEM model of the glass fiber and the epoxy matrix are carried out to validate the main mechanical properties, as listed in Table .1. Figures 7 and8 show the respective strain-stress relationships, where both the longitudinal and transverse strains are recorded against the applied tensile load. The parameters in the above calibrated DEM models are then incorporated in the laminae model satued in the following sections. In this study transverse cracking of a glass/epoxy lamina under tensile loading, as shown in Fig. 9, is studied, where three types of fiber distribution and three fiber volume fractions are considered. The interface strength is smaller than that of the matrix, with both the tensile and shear strength being 50 MPa, and the interface fracture energy is chosen as 100 J/m 2 . A large length-to-depth ratio is used to eliminate the far field loading effect. The particles with a minimum radius of R min= 5 µm (R max =1.5R min ) are generated first within a rectangular square of 1.5 mm × 1.0 mm in the middle of the specimen, where the fibers are inserted, to represent the composite. The central part is then extended to the both sides by a sequence of four rectangular assemblies composed of larger matrix particles. The final size of the model is 1.5 mm × 5 mm. The fibers are located in a small square area of 650 µm × 650 µm in the middle of the specimen in order to create a fiber reinforced composite laminae. An initial crack of length a was created in the matrix at the top of the mid-span by deleting the corresponding particle elements. This is used to ensure that a transverse crack would propagate towards the fiber reinforced area to simulate the experimental study in [START_REF] Vejen | Transverse crack growth in glass/epoxy composites with exactly positioned long fibres. Part I: experimental[END_REF]. A constant displacement loading with a velocity of 0.1m/s is applied at the both ends of the specimen.

In order to study the effect of fiber distribution on the transverse cracking, rectangular, hexagonal and random distributions of 49 fibers are used in the DEM models as shown in Fig. 10. The numerical results are compared with the experimental observations of transverse cracking path in Fig. 11. As the load applied at the two ends increases, matrix cracking, which is initiated from the preset crack, propagates into the interior area of the material and subsequently interacts with interfacial debonding. The cracks finally propagate through the entire matrix as illustrated in Fig. 10 (i) and Fig. 10 (ii). In all of the three cases of fiber distribution, when the force between two matrix particles exceeds the bond strength of the material, the matrix exhibits plastic deformation (represented by the black lines in Fig. 10) and finally breaks when the fracture energy is reached. At this moment, the black lines are replaced by the red lines in Fig. 10 to show the formation of cracks. Similarly, fiber/matrix interface started to yield when its strength is exceeded (represented by the brown lines in Fig. 10), and debonding occurs when its strength is softened to zero, represented then also by the red lines. From Fig. 10, it can be seen that the entire dynamic damage process of the initiation and propagation of cracking can be naturally simulated by the DEM model.

Results and Discussions

It can be found from the DEM results that for all the three fiber distributions the transverse cracks propagate through the middle part of the fiber reinforced area, and the crack development is altered and/or interrupted by fibers and fiber/matrix interfaces, resulting in a zigzag path, as shown in Fig. 10 (ii). Plastic deformation of the matrix and fiber/matrix interfacial yielding can be observed on both sides of the transverse cracks, but more severe plastic deformation of the matrix and fiber/matrix interfacial yielding are found in the vicinity of the transverse crack. In the laminae with rectangular distribution of fibers, interfacial debonding always occurs when the crack propagates toward the fiber/matrix interface, though some of the debondings occurs on the left hand side of a fiber, while others takes place on the right hand side. A similar phenomenon has been recorded in the experimental investigation [35], as shown in Fig. 11(a).

In the laminae with hexagonal arrangement of fibers, the transverse cracking path is slightly different. It joins the fiber/matrix interfacial debondings as well as penetrates through the matrix between two horizontally neighboring fibers, and directly propagates forward, as highlighted in zone (A) of The propagation of transverse cracking in randomly distributed fiber reinforced laminae is much more complex in reality, and is successfully simulated by the DEM model without introducing any assumptions on how the crack should propagate, as shown by the numerical results in Fig. 10c. The direction of the transverse cracking is more or less random but generally followed the middle line of the specimen. It is also a mixture of joining interfacial debonding and matrix cracking.

It is evident that the DEM results have provided a unique insight of the microscopic damage within the composite laminae. One of the unique features of the DEM model is that the dynamic simulation can record all the damage accumulated during the loading process. From these records, it can be found that the distribution of fibers not only affects the transverse cracking path, but also causes different residual damage, such as matrix plastic deformation and fiber/matrix interfacial yielding in the material. Fig. 12 presents a statistic summary of the damage accumulated in the loading process. In the laminae with rectangular fiber distributions, the number of residual damage, in terms of the total number of the particle contacts exhibiting either matrix plastic deformation or fiber/matrix interface yielding, is smaller than that with hexagonal or randomly distributed fibers. This suggests that laminae with rectangular distribution of fibers can offer better residual mechanical properties. This is due to the fact that a rectangular distribution of fibers can offer more uniform stress distribution within the laminae, resulting in a lower level of stress concentration, and then a reduction of plastic deformation and fiber/matrix interface debonding. For randomly distributed fiber reinforcement, Figures 13 and14 show the results obtained by the DEM model with different fiber volume fractions. In the laminae with lower fiber volume fractions, a more diverted transverse cracking path is observed, probably due to the fact that more matrix material between two neighboring fibers has left more freedom for the transverse crack to pass through, which creates a relatively less uniform path. When the fiber volume fraction is higher, the propagation of the transverse cracking is largely guided by fiber/matrix interfacial debonding, and therefore, the path of transverse cracking is more direct, resulting in a much smoother crack.

Also in the case of random distribution of fibers, it can be found that the total fiber/matrix interfacial yielding increases, while the number of matrix plastic deformation decreases, with the increase of the fiber volume fraction, as illustrated in Fig. 14. This is due to the reduction of matrix volume and average distance between two neighboring fibers, making the interfaces much easier to yield.

Conclusions

Discrete element method (DEM) has been employed to simulate transverse cracking in composite laminae, due to its intrinsic advantages in modeling microscopic damage and fracturing in multiphased materials. Three types of fiber distributions, i.e., rectangular, hexagonal and random distributions, have been simulated to study the effect of fiber distribution on the transverse cracking. The initiation as well as dynamic propagation of transverse cracking has been captured by the DEM model, showing good agreement to the experimental observations. Furthermore, the DEM simulations have provided unique insights of the microscopic cracking and damage in forms of matrix plastic deformation and fiber/matrix interface yielding. The effect of fiber volume fraction was also studied by using different fractions in the DEM modeling for randomly distributed fibers.

It is a unique feature of the DEM that allows a quantitative study of the residual microscopic damage within a fractured material. This has led to the findings that the distribution and volume fraction of fibers not only affect the transverse cracking path, but also the mechanical behaviors of the material through the populations of the residual matrix plastic deformation and fiber/matrix interfacial yielding.

It can be concluded that the DEM developed in this paper is particularly suitable for the modeling of complex multi scale damage and cracking in multiphase composite materials, and can be further modified to include thermal effect and applied to model more complex composite systems in the future, such as cross-ply laminates.

Acknowledgement

The authors wish to thank Professor N L McCartney, the National Physical Laboratory, for reviewing an early version of the manuscript and providing detailed criticism and comments. The second author would like to acknowledge the ORSAS scholarship and the University of Leeds studentship for the financial support of his PhD study. The third author would like to thank the NSFC of China for the financial support of this collaboration research under the grant of No.50875224.

Fig. 1 Fig. 3

 13 Fig.1 Random arrangement of particles

Fig. 4

 4 Fig.4 Convergence of loading velocities for quasi-static simulation

Fig. 5 Fig. 6

 56 Fig.5 DEM model of DCB test

Fig. 7 Fig. 8

 78 Fig.7 Stress-strain curve in DEM simulation of fiber tensile test

Fig. 10 DEM

 10 Fig.10 DEM dynamic simulation of transverse cracking in laminae ('a', 'b' and 'c' represent rectangular, hexagonal and random distribution of fibers, respectively; Red, black and brown lines describe transverse cracking, matrix plastic deformation and interface yield, respectively.)

Fig. 10 bFig. 11

 1011 Fig.11 Experimental observation of cracks in laminae with different arrangement of fibers in [35]: (a) rectangular; (b) hexagonal

Fig. 12

 12 Fig.12 Residual damage in laminae with different fiber distributions

Fig. 13

 13 Fig.13 Transverse cracking in laminae with random fiber distribution and different fiber volume fractions: (a) 23.2%, (b) 37.6% and (c) 56.2%.

Fig. 14

 14 Fig.14 Residual damage in laminae with different fiber volume fractions

Table 1

 1 Material properties of the DEM model

	Fiber modulus (glass)

σε x σε y

(a) (b) (c)


Nelma Gomes
email: nelmagomes@deb.uminho.pt

Yves Waché
email: yves.wache@u-bourgogne.fr

José A Teixeira
email: jateixeira@deb.uminho.pt

Isabel Belo
email: ibelo@deb.uminho.pt

Microbial and Enzyme Technology

Keywords: castor oil, -decalactone, droplet size distribution, emulsions, methyl ricinoleate

Oil-in-water emulsions composed of methyl ricinoleate (MR) or castor oil (CO) as the organic phase, stabilized by Tween 80, are in the basis of the biotechnological production of -decalactone. Yarrowia lipolytica was used due to its ability to grow on hydrophobic substrates and to carry out the biotransformation. The characterization of oil droplets size distribution by laser granulometry was performed under different oil concentrations. The impact of the presence of cells on droplets size was also analyzed as well as the relevance of washing inoculum cells. Furthermore, the granulometric characterization of the emulsions was related with -decalactone production and it was observed that, in the presence of non-washed cells, the smaller droplets disappeared, using both oils, which increaseddecalactone concentration. This suggests that the access of cells to the substrate occurs by their adhesion around larger oil droplets.

Introduction

The biotechnological production of -decalactone, a fruity-like aroma compound, results from the biotransformation of ricinoleic acid and it can be carried out by several microorganisms, among which the dimorphic yeast Yarrowia lipolytica (Aguedo et al. 2005a).

The medium usually used is composed by an oil-in-water emulsion stabilized by a non-ionic surfactant, Tween 80. The oil is the source of ricinoleic acid, usually methyl ricinoleate (MR) or castor oil (CO).

The growth of cells in a medium containing lipids occurs directly on droplets as their solubility in water is limited [START_REF] Gutierrez | Hydrocarbon uptake in hydrocarbon fermentations[END_REF] and it is, thus, influenced by the size of the oil droplets [START_REF] Bakhuis | Correlation between growth of Candida lipolytica and size of droplets in the hydrocarbon containing medium[END_REF].

The interfacial surface of the oil droplets and consequently the surface between the two liquid phases is a determining factor for the degradation of the hydrophobic substrate and, thus, for the cells, growth and aroma production (Aguedo et al. 2005b). The lipidic droplets size is influenced by numerous parameters which may be related to physical and chemical characteristics of the medium (pH, ionic strength, presence of surfactants, etc) [START_REF] Wilde | Interfaces: their role in foam and emulsion behavior[END_REF] or to the presence of microorganisms (inoculum concentration, cells surface properties, etc). During the biotransformation of ricinoleic acid intodecalactone, there is a direct contact between the cells surface and the substrate droplets which may occur by adhesion of cells to the larger droplets surface or by adsorption of micro-droplets at the cells surface [START_REF] Bakhuis | Correlation between growth of Candida lipolytica and size of droplets in the hydrocarbon containing medium[END_REF].

Since the access of the cells to the substrate depends on the emulsion characteristics, emulsions were characterized by determining the oil droplets size distribution, through the laser granulometry technique, using both oils (MR and CO), at different concentrations and in the presence or absence of cells.

Moreover, since the cells surface properties may influence the interaction between the oil and cells, different methodologies were used for inoculum cells preparation, according with [START_REF] Gomes | The use of methyl ricinoleate in lactone production by Yarrowia lipolytica: Aspects of bioprocess operation that influence the overall performance[END_REF].

This procedure had an impact on cells hydrophobic character and the effect on emulsion behaviour was analyzed in the present work. In a last instance, the impact of each emulsion on -decalactone production was also investigated.

Materials and methods

Microorganism, media and culture conditions

Yarrowia lipolytica W29 (ATCC20460) was cultured for 48 h on YPDA medium (30 g agar l -1 , 20 g glucose l -1 , 20 g peptone l -1 , 10 g yeast extract l -1 ) at 27 ° C and used to inoculate (to reach OD 600 of ~0.25) a 500 ml baffled Erlenmeyer flask containing 200 ml glucose medium (YPD medium: glucose 20 g glucose l -1 , 20 g peptone l -1 , 10 g yeast extract l -1 ). Flasks were shaken at 140 rpm and 27 °C for 19 h until the cultures reached the late growth phase. This suspension was used to inoculate the biotransformation medium.

Two sources of ricinoleic acid were used: methyl ricinoleate (MR) and castor oil (CO) at different concentrations. The composition of the biotransformation medium was: 6.7 g YNB (yeast nitrogen base) with amino acids l -1 , 2.5 g NH 4 Cl l -1 , from 10 to 50 g MR or CO l -1 and from 1 to 5 g Tween 80 l -1 . Two methodologies were tested for inoculation of the biotransformation medium: with or without the cell washing step. In the first case, cells from the YPD medium were washed three times in sterile water, centrifuged (6000 g, 5 min) and transferred to the biotransformation medium. In the second, the components of the biotransformation medium were directly added to the YPD medium containing the cells [START_REF] Gomes | The use of methyl ricinoleate in lactone production by Yarrowia lipolytica: Aspects of bioprocess operation that influence the overall performance[END_REF].

All chemicals were purchased from Sigma-Aldrich except for CO, which was purchased from José M.

Vaz Pereira (Sintra, Portugal) and MR, which was kindly supplied by Stéarinerie Dubois (Boulogne, France).

-Decalactone extraction and quantification -Decalactone was extracted from 2 ml broth with 2 ml diethyl ether. The organic phase was analyzed by GC as described in previous work [START_REF] Gomes | The use of methyl ricinoleate in lactone production by Yarrowia lipolytica: Aspects of bioprocess operation that influence the overall performance[END_REF].

Laser granulometry

The droplet size distribution of MR and CO emulsions was evaluated using a Malvern Mastersizer laser granulometer (Model S2-01, Malvern Instruments, Malvern, UK), based on the Fraunhofer approximation [START_REF] Agrawal | Laser diffraction size analysis[END_REF].

Emulsions with MR and CO were analyzed with and without cells precultured during 19 h on YPD medium. For the essays without cells, the emulsions were shaken during 1.25 h prior to analysis. For the experiments with cells, the emulsions were shaken for 1 h before inoculation and then for more 15 min, with cells, before analysis. A cellular suspension prepared in sterile water was also analyzed, in order to determine the size of the yeast cells. All assays were performed in triplicate.

Results and discussion

The granulometric analyses of the emulsified MR and CO media, without cells, indicated the existence of two distinct droplet populations (Figure 1). In the population with smaller dimension, the droplets size remained practically with the same mean value among all MR (1.5 m) and CO (1.9 m) concentrations tested. In the larger size population, the mean droplets diameter increased from 61 m to 96 m, with the increase in MR concentration, and from 27 m to 64 m, with the increase in CO concentration.

Moreover, there was an increase in percentage of larger droplets with the oil concentration increase.

(Figure 1) [START_REF] Aguedo | Increased electron donor and electron acceptor characters enhance the adhesion between oil droplets and cells of Yarrowia lipolytica as evaluated by a new cytometric assay[END_REF] have also characterized an emulsion composed by 10 g MR l -1 . The results obtained by these authors differ from those presented in this work since they have reported the existence of three distinct droplet populations with mean sizes of 0.1 m, 3 m and 43 m, in the MR emulsion.

These divergences may be justified by a variation in the concentration of Tween 80 (0.2 g l -1 vs. 1 g l -1 , in their work and in the present one, respectively). Thus, it seems that the use of a higher concentration of surfactant has an impact upon the MR droplets size, since the smaller droplets (0.1 m) disappeared and a population of larger droplets emerged: 61 m in this study, compared with 43 m in their work.

Prior to the determination of droplets size in the presence of cells, the cells size distribution was evaluated by suspending the cells in sterile water. A population with mean values of 5.4 µm was detected, corresponding to an oval morphology [START_REF] Lopes | Yarrowia lipolytica lipase production enhanced by increased air pressure[END_REF]).

(Figure 2) When cells were present at around 1.5 10 9 ml -1 in the MR and CO emulsions, there were still two populations (in some cases, three) with different particle sizes (Figures 3 and4). However, drops with 1.5-1.9 µm size disappeared in all emulsions. This disappearance, due to the presence of cells, may be explained by the adhesion of smaller substrate droplets to the cells surface and to the higher number of cell particles that would overshadow the smallest oil droplets [START_REF] Aguedo | Increased electron donor and electron acceptor characters enhance the adhesion between oil droplets and cells of Yarrowia lipolytica as evaluated by a new cytometric assay[END_REF]). Nevertheless, the population of cells remained with a mean size of approximately 5.4 µm. The existence of this population was also reported by [START_REF] Aguedo | Increased electron donor and electron acceptor characters enhance the adhesion between oil droplets and cells of Yarrowia lipolytica as evaluated by a new cytometric assay[END_REF].

(Figures 3 and4) When washed cells were used (Figure 3), a population had a mean oil droplet size of 0.9 µm, although in a small percentage number, in all experiments and for both oils. This decrease in the oil droplets mean size with the addition of Y. lipolytica washed cells was also stated by [START_REF] Waché | Medium-size droplets of methyl ricinoleate are reduced by cell-surface activity in the -decalactone production by Yarrowia lipolytica[END_REF], when using emulsions composed by 5 g MR l -1 and 1 g Tween 80 l -1 .

Another population, with a mean droplet size comprised between 32 µm and 40 µm, was also noticed in all experiments and for both substrates, although with a larger number of particles in CO emulsions. [START_REF] Aguedo | Increased electron donor and electron acceptor characters enhance the adhesion between oil droplets and cells of Yarrowia lipolytica as evaluated by a new cytometric assay[END_REF] also described the existence of a population of MR droplets with a mean size of 42 µm, in emulsions with Y. lipolytica washed cells.

In the experiments involving non-washed cells (Figure 4), as the concentration of MR increased, the size of oil droplets also increased slightly, ranging from 40 µm to 69 µm, except for the experiment with 50 g MR l -1 , which presented a mean size of 52 µm.

When comparing the experiments with washed and non-washed cells, using both oils, one can observe that the smallest oil droplets disappeared (except for the case with 50 g CO l -1 , where the existence of a small population of oil droplets with a mean size of 0.8 µm was noticed). This disappearance is probably related to the higher hydrophobicity of non-washed cells, as proved in previous work [START_REF] Gomes | The use of methyl ricinoleate in lactone production by Yarrowia lipolytica: Aspects of bioprocess operation that influence the overall performance[END_REF], allowing droplets to adsorb to their surface, turning its presence not measurable.

Table 1 shows that the accumulation of -decalactone is directly dependent on the MR concentration, with the greater production of aroma occurring when the substrate concentration is higher. The accumulation of -decalactone was also greater when non-washed cells are used. In this situation, the MR droplets size is higher (the smaller population with 0.9 µm disappeared), suggesting that higher particle sizes are more favorable to the production of the aroma, indicating that the adhesion of cells to the substrate occurs by the linkage of cells around larger oil droplets [START_REF] Bakhuis | Correlation between growth of Candida lipolytica and size of droplets in the hydrocarbon containing medium[END_REF]. In fact, cells possess a much smaller size than the oil droplets and the number of cell particles is far superior to the number of oil particles, which seems to corroborate the hypothesis proposed.

(Table 1) With respect to the use of CO and non-washed cells, the size of oil droplets varied from 34 µm to 46 µm, for the larger population (Figure 4).

When castor oil is used as substrate of the process, there is no aroma production when washed cells are used in the biotransformation. When using non-washed cells, the most adequate oil concentration to be employed, in terms of overall productivity, is 30 g l -1 , suggesting that the most favorable droplets size corresponds to a diameter of approx. 46 µm, the larger size determined among all CO emulsions. This corroborates the hypothesis previously suggested about larger particle sizes being more favorable to the aroma production.

The main differences obtained in the aroma production under different operating conditions cannot be attributed solely to the size of oil droplets; the quantity of available substrate for cells to use in the biotransformation, as well as to the cells surface properties, which determines how the substrate is assimilated, need to be taken into account. [START_REF] Gomes | The use of methyl ricinoleate in lactone production by Yarrowia lipolytica: Aspects of bioprocess operation that influence the overall performance[END_REF] verified, by means of the MATH test, that there was a greater adhesion of non-washed cells to hexadecane, than of washed cells. These differences point out that non-washed cells are more compatible with the substrate, which might also explain the divergences obtained in the production of -decalactone, apart from the differences in the oil droplets size.

Conclusion

This study aimed to characterize the oil droplets size distribution on different oil-in-water emulsions, composed by different concentrations of methyl ricinoleate or castor oil and stabilized by a non-ionic surfactant. The impact of the presence of cells on droplets size was shown, as well as the relevance of a washing step of the inoculum cells.

In the experiments without cells, an increase in both oil concentrations leaded to an enlargement of the mean oil droplets size. In the presence of cells, there are two or three populations with different particle sizes. However, when cells were washed, a smaller population of both oil droplets was detected.

Moreover, the granulometric characterization of the emulsions was related with -decalactone accumulation profiles and allowed to conclude that larger oil droplets favors the aroma production, suggesting that, in this case, the access of cells to the substrate occurs by their adhesion around larger oil droplets.

The obtained results are an important contribution for a better understanding of the factors affecting the assimilation of hydrophobic compounds by Yarrowia lipolytica, clearly demonstrating the impact of substrate concentration, as well as of cells, on droplet size distribution and, consequently, ondecalactone production.

Figure Captions

Fig. 1 Fig. 2 Fig. 3 Fig. 4

 1234 Fig.1Methyl ricinoleate (MR) and castor oil (CO) droplets size distribution (in µm) in emulsions without cells, related to the number of particles (in %). Data are the mean of three independent experiments

Acknowledgments

8

The authors acknowledge Fundação para a Ciência e Tecnologia (FCT) for the financial support provided (SFRH/BD/28039/2006) and CRUP (Cooperation project F-41/06).

Fig. 2


Suman Kundu

Yurii S Aulchenko

Cornelia M Van Duijn

A Cecile

J W Janssens
email: a.janssens@erasmusmc.nl

PredictABEL: an R package for the assessment of risk prediction models

Keywords: risk prediction, genetic, assessment, measures, software

The rapid identification of genetic markers for multifactorial diseases from genomewide association studies is fuelling interest in investigating the predictive ability and health care utility of genetic risk models. Various measures are available for the assessment of risk prediction models, each addressing a different aspect of performance and utility. We developed PredictABEL, a package in R that covers descriptive tables, measures and figures that are used in the analysis of risk prediction studies such as measures of model fit, predictive ability and clinical utility , and risk distributions, calibration plot and the receiver operating characteristic (ROC) plot.

Tables and figures are saved as separate files in a user-specified format, which include publication-quality EPS and TIFF formats. All figures are available in a ready-made layout, but they can be customized to the preferences of the user. The package has been developed for the analysis of genetic risk prediction studies, but can also be used for studies that only include non-genetic risk factors. PredictABEL is freely available at the websites of GenABEL (http://www.genabel.

Introduction

The rapid identification of genetic markers for multifactorial diseases from genomewide association studies is fuelling interest in investigating the predictive ability and health care utility of genetic risk models. Genetic risk models are investigated for their potential to target diagnostic, preventive and therapeutic interventions for multifactorial diseases. Implementation of these models in health care requires a series of studies that encompass all phases of translational research [START_REF] Khoury | The continuum of translation research in genomic medicine: how can we accelerate the appropriate integration of human genome discoveries into health care and disease prevention?[END_REF][START_REF] Hlatky | Criteria for evaluation of novel markers of cardiovascular risk: a scientific statement from the American Heart Association[END_REF], starting with a comprehensive evaluation of genetic risk prediction.

Various measures are available for the assessment of risk prediction models, each addressing a different aspect of performance and utility [START_REF] Mcgeechan | Assessing new biomarkers and predictive models for use in clinical practice: a clinician's guide[END_REF][START_REF] Steyerberg | Assessing the performance of prediction models: a framework for traditional and novel measures[END_REF]. The GRIPS Statement recommends that transparent and complete reporting should provide a description of the risk factors and the risk model by reporting univariate and multivariate odds ratios for the predictors, present risk distributions for individuals with and without the outcome of interest, and report measures of model fit, predictive ability and others, if pertinent [START_REF] Janssens | Strengthening the reporting of genetic risk prediction studies: The GRIPS Statement[END_REF][START_REF] Janssens | Strengthening the reporting of genetic risk prediction studies (GRIPS): elaboration and explanation[END_REF]. Examples of measures include the Hosmer-Lemeshow statistic [START_REF] Hosmer | A comparison of goodness-of-fit tests for the logistic regression model[END_REF] and Nagelkerke's R 2 [START_REF] Nagelkerke | A note on a general definition of the coefficient of determination[END_REF] for model fit, the area under the receiver operating characteristic (ROC) curve (AUC) [START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (ROC) curve[END_REF] and integrated discrimination improvement (IDI) [START_REF] Pencina | Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond[END_REF] for predictive ability, and percentages of total reclassification [START_REF] Cook | Use and misuse of the receiver operating characteristic curve in risk prediction[END_REF] and net reclassification improvement (NRI) [START_REF] Pencina | Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond[END_REF] for clinical utility.

Even though the assessment of risk prediction models is relatively standard, there is no single statistical package that would allow for the computation and production of all these measures and plots. Therefore, we developed PredictABEL, a freely available R package, which contains functions to obtain all descriptive tables, measures and plots that are used in genetic risk prediction studies.

Description of PredictABEL

The core part of PredictABEL comprises functions for the assessment of risk prediction models. The measures and plots covered in PredictABEL are listed in Table 1. Most functions can be applied to predicted risks, risk scores or any other continuous predictor variable, but some to predicted risks (probabilities) only.

Predicted risks and genetic risk scores can be obtained using functions in the package, but they can be imported from other programs as well. The functions to obtain predicted risks using logistic regression analysis are specifically written for models that include genetic variables, eventually in addition to non-genetic factors, but they can also be applied to construct models based on non-genetic risk factors only.

Genetic risk scores can be computed as unweighted and weighted risk scores, where weights are obtained from uploaded data or imported from meta-analyses, e.g., as beta coeffcients.

The tables and plots generated using PredictABEL are saved as separate files in the working directory. Tables can be saved as Excel or tab-delimited text files and figures can be saved as publication-quality EPS or TIFF files or as JPEG files for insertion in manuscripts. All figures are available in a ready-made layout, but they can be customized to the journal style or preferences of the user. A hypothetical dataset and examples of use are included in the package to demonstrate all functions.

Example

The hypothetical dataset included in the package was reconstructed from an empirical study on age-related macular degeneration (AMD) [START_REF] Seddon | Prediction model for prevalence and incidence of advanced age-related macular degeneration based on genetic, demographic, and environmental variables[END_REF], using a simulation method that has been described in detail elsewhere [START_REF] Janssens | Predictive testing for complex diseases using multiple genes: fact or fiction?[END_REF]. Based on published frequencies and odds ratios of the genetic variants and non-genetic risk factors implicated in AMD and on published population disease risks, we created a dataset that contains genotype data and disease status for 10,000 individuals. Predicted risks were obtained using logistic regression analysis, for which the codes are provided in the package. Two risk models were constructed: a model based on non-genetic risk factors only and a model based on genetic and non-genetic predictors.

Figure 1 presents three examples of plots that are produced by PredictABEL.

Figure 1a shows distributions of predicted risks based on genetic and non-genetic factors for individuals with and without AMD. The degree of overlap between the two histograms is indicative for the discriminative accuracy of the risk model. This discriminative accuracy is assessed by the AUC and visualized in a ROC plot. Figure 1b presents the ROC curves for the two risk models. The figure shows that the model with genetic factors had a higher AUC than the model without. Using the same function, the AUC values were quantified as 0.80 and 0.74. Finally, Figure 1c presents the calibration plot for the risk model based on the genetic and non-genetic variables as predictors, which shows how well predicted risks match observed risks.

The calibration plot suggests that the model was well calibrated, which was supported by the non-significance of the Hosmer-Lemeshow test (p = 0.65).

Finally, Table 2 presents an example of the reclassification table and statistics that are produced by PredictABEL. The reclassification table presents the categorization into risk groups according to the initial and updated risk models. The table provides information about the total number of individuals that change between risk categories and about correct and incorrect reclassification. The percentage of total reclassification and NRI are calculated from the reclassification table. The table indicates that net 8.8% of the individuals without AMD and 9.6% of those with AMD would be correctly reclassified when the clinical model was updated by the addition of genetic factors.

Conclusions

PredictABEL is a comprehensive software package, designed for the development and assessment of genetic risk prediction models. PredictABEL is a part of the GenABEL software suite for statistical genomics [START_REF] Aulchenko | GenABEL: an R library for genome-wide association analysis[END_REF][START_REF] Aulchenko | ProbABEL package for genome-wide association analysis of imputed data[END_REF] and for that reason written in R to enable easy transfer of data from gene discovery to genetic prediction studies. A detailed manual is available that demonstrates and explains all the functions in the package. The manual is accessible for researchers who do not regularly use R software. The manual and the package are freely available from the GenABEL project website (http://www.genabel.org) and from CRAN (http://cran.r-project.org/).

The current version of PredictABEL (version 1.1) includes all basic descriptive tables, measures and plots that are used in the assessment of risk prediction models. Planned extensions of the package include other strategies to construct risk models, e.g., using Cox Proportional Hazards analysis for prospective data, and functions to construct simulated data for the evaluation of genetic risk models [START_REF] Janssens | Predictive testing for complex diseases using multiple genes: fact or fiction?[END_REF]. Furthermore, we will optimize the interconnectivity between PredictABEL and other packages in the GenABEL suite.

Where the GRIPS Statement aims to improve the transparency, quality and completeness of reporting [START_REF] Janssens | Strengthening the reporting of genetic risk prediction studies: The GRIPS Statement[END_REF][START_REF] Janssens | Strengthening the reporting of genetic risk prediction studies (GRIPS): elaboration and explanation[END_REF], PredictABEL has similar goals for the assessment of genetic risk prediction studies. The collection of all measures and plots in a single, software package gives a comprehensive overview of the various measures that are available for the assessment of risk prediction studies. This overview emphasizes that different measures are available to answer different questions in the assessment of risk models and facilitates the selection of the most appropriate measure for the question under study.

Figure 1

 1 Figure 1 Example graphs produced by PredictABEL.

Acknowledgements

This work was supported by the Vidi grant from the Netherlands Organization for Scientific Research (NWO), the Young Investigator grant from the Erasmus University Medical Center Rotterdam and by the Center for Medical Systems Biology within the framework of the Netherlands Genomics Initiative.

* These functions can only be used when the logistic regression model is constructed using the functions in PredictABEL.