Carlos Murillo

Ghislain Irakoze

Karine De Oliveira Vigier

Michel Delmas

François Jérôme

Yolande Peres-Lucchese

Martine Urrutigoity

Patrick Cognet
email: patrick.cognet@ensiacet.fr

Yolande Pérès

Modeling of Ethylene Glycol Production from Glucose in a Semi-Continuous Reactor

Keywords: Ethylene glycol, glucose, reaction kinetics, modelling, semi-continuous reaction

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Ethylene glycol is a commodity whose consumption is associated with various industries such as energy, chemicals, automotive, textiles, transportation, and manufacturing technologies [1]. For instance, it can be used as an antifreeze in automobile radiators due to its ability to lower the freezing point when mixed with water. Additionally, it is widely utilized as a raw material for the manufacture of polyester fibers [START_REF] Rebsdat | Ullmann's Encyclopedia of Industrial Chemistry[END_REF]. For this reason, the global production and consumption of this diol were about 20 million metric tons in 2010 and reached 25 million tons in 2015 [1], [START_REF] Zhao | [END_REF]. Several synthesis routes were developed to produce EG [1], [START_REF] Rebsdat | Ullmann's Encyclopedia of Industrial Chemistry[END_REF], but the major part of the market share relies on the production from ethylene oxide hydration. This fact associates the production of this diol to the petroleum industry, which represents an economic manufacture technology. the diminution of the fossil-fuel resources and the increasing environmental concerns make it necessary to consider other economical and energy-efficient processes [1]. At this moment, the manufacture of ethylene glycol from other sources such as syngas [4] or biomass [5] become more interesting alternatives.

Biomass is a convenient substitute for fossil resources because it can be transformed into useful organic compounds [6]. Indeed, glucose can be considered as a potentially promising feedstock for the production of EG because it is available for continuous processing on a large scale [7], [8]. For this purpose, agricultural wastes can be treated to obtain cellulose, which can subsequently be hydrolyzed to glucose. Then, glucose is transformed into glycolaldehyde with a homogeneous catalyst. Finally, the aldehyde is hydrogenated to ethylene glycol with a solid catalyst dispersed in the mixture.

The conversion of glucose into ethylene glycol constitutes an opportunity to benefit from a renewable raw material and generate a massive reduction in CO2 emission [9]. However, the main drawback of this synthesis route is associated with the presence of several side reactions that decrease considerably the yield of the chemical process. For this reason, this technique requires a thorough analysis for determination of the operating conditions that will enhance the conversion of glucose into ethylene glycol. Previously, Zhao et al. developed a kinetic model for the main reactions of this process with an approach based on power law equations with apparent constants [START_REF] Zhao | [END_REF]. In this manner, their mathematical analysis evaluated the influence of the reaction temperature and the feed rate on the yields obtained in a semi-batch reactor. These results allowed concluding that the selectivity to ethylene glycol is increased when the reaction temperature is between 200 and 240 °C. Additionally, Zhang et al. characterized the kinetics of the glucose and glycolaldehyde hydrogenations occurred at the surface of the solid catalyst by establishing a set of Langmuir-Hinshelwood-Hougen-Watson (LHHW) equations for temperatures between 100 and 130 °C [8].

The aim of this study is to propose a global kinetic model based on the combination of the approaches based on power law models for the homogeneous reactions and LHHW models for the heterogeneous reactions. For this purpose, the parameters of the LHHW equations are estimated according to literature data for the temperature range suggested by Zhao et al. [START_REF] Zhao | [END_REF], [5]. The new model is then implemented in a semi-batch reactor mass balance. In this manner, the influence of operating conditions of a semi-batch operation such as temperature, hydrogen partial pressure, catalyst concentrations and glucose feeding rate, on EG yield and selectivity can be described. Subsequently, a sensitivity analysis of the key operating variables is performed. These results complement the previous experimental and computational studies of the glucose conversion into ethylene glycol within a semi-continuous reactor.

Chemical reaction network

The high chemical activity of glucose allows performing its catalytic conversion to ethylene glycol at mild conditions with high efficiencies [5]. In addition, this compound can be obtained from highconcentration sugar solutions [START_REF] Schenck | Ullmann's Encyclopedia of Industrial Chemistry[END_REF]. Therefore, ethylene glycol is produced from this monosaccharide at high concentrations and low energy consumption in the ethylene glycol purification. However, the synthesis route also has certain drawbacks due to the occurrence of competitive reactions that promote the conversion of glucose to hexitols, polyols with 3 or 4 carbon atoms, humins, and gases. Hence, it is necessary to describe the ethylene glycol production as a part of a reaction network in order to identify the most appropriate operating conditions. The six types of chemical reactions that compose the chemical process of this study are shown in in supporting information and described as follows:

i. Retro-aldol condensation: This chemical reaction corresponds to the decomposition of a β-hydroxyl carbonyl into an aldehyde or ketone. It is observed in this chemical process following the conversion of glucose into erythrose and glycolaldehyde. Additionally, it is associated with the formation of dihydroxyacetone (DHA) and 2,3 dihydroxypropanal (DHP) from fructose [START_REF] Zhao | [END_REF], [START_REF] Zhang | [END_REF]. Retro-aldol condensation is affected by the medium conditions. The reaction can be carried out in an acid or basic medium and in supercritical water [12]. According to literature, several tungstic compounds such as tungsten oxide, tungstophosphoric acid, tungsten acid, tungstosilicic acid and ammonium metatungstate (AMT) can be used for this purpose [5], [13], [14]. AMT is usually selected since it provides the highest ethylene glycol yields due to its solubility in aqueous media [15]. According to this analysis, the kinetic model developed in this study considers the addition of AMT as the homogeneous catalyst. The use of AMT reduces both glucose and glycolaldehyde hydrogenation reaction rates [8]. Hence, the global ethylene glycol yield is notably enhanced by the addition of this homogeneous catalyst in spite of the depletion of the glycolaldehyde hydrogenation.

ii.

Hydrogenation: This reaction is mainly associated with the formation of ethylene glycol from glycolaldehyde. However, the hydrogen dissolved in the aqueous solution also reacts with saccharides (glucose, fructose, and mannose) to produce hexitols (sorbitol and mannitol). Furthermore, erythrose, DHA and DHP can also produce three-carbon or four-carbon polyols due to competitive hydrogenation. These compounds are highly stable; therefore, they cause the diminution of the global ethylene glycol yield [START_REF] Zheng | Reaction Pathways and Mechanisms in Thermocatalytic Biomass Conversion I: Cellulose Structure, Depolymerization and Conversion by Heterogeneous Catalysts[END_REF].

Hydrogenation reactions are carried out with the dispersion of a heterogeneous catalyst. Nickel compounds are commonly used for industrial-scale processes due to their low costs. However, due to their several drawbacks [START_REF] Aho | [END_REF], these catalysts are currently being replaced by others whose structures are composed of a noble metal supported on activated carbon (AC). Zhao et al. developed a comparative analysis of the yields obtained with various catalysts: Ru/AC, Pt/AC, Pd/AC, Ir/AC and Raney Ni [5]. Their experimental tests evidenced a higher yield and selectivity to ethylene glycol with ruthenium (4% Ru/AC). Besides this advantage, this metal-based catalyst is characterized by a comparatively low price with regard to other noble metals [9]. This study is performed with experimental data obtained in the literature for 4% Ru/AC as heterogeneous catalyst [START_REF] Zhao | [END_REF], [5].

iii.

Thermal side reactions: Some insoluble by-products (humins or tars) can be produced due to high reactant concentrations and reaction temperatures over 260 °C. Thus, it is necessary to control the temperature below this value.

iv.

Isomerization: Glucose injected into the aqueous solution can isomerize into fructose. This reaction is enhanced by the presence of basic active carbon [14]. Fructose will produce different by-products by subsequent hydrogenation reactions. In fact, fructose produces DHA and DHP which can transform respectively into 1,2 propylene glycol and glycerol.

v.

Epimerization: Glucose can be converted into mannose due to a rotation of the C2-C3 bond. Subsequently, mannose can be directly hydrogenated to mannitol (hexitol). The epimerization reaction requires a substantial reorganization of the water shell of the molecule in the aqueous solution. For this reason, the mannose production rate is slower than that of fructose [18]. However, this reaction rate can be increased by the presence of basic catalysts and Lewis acids [19].

vi.

Hydrogenolysis: Ethylene glycol can decompose due to a reaction with dissolved hydrogen to produce methane [20]. This transformation represents another issue of an overheating process in the chemical reactor.

After describing the different reactions involved in the glucose conversion into ethylene glycol, a kinetic model is developed in the next section, based on the simplification of the reaction network.

3

Kinetic model

The influence of the catalyst concentrations and the hydrogen partial pressure on the reaction yields can be analyzed numerically with the development of a kinetic model that takes into account these operating parameters. For this purpose, some simplifications are proposed for the reaction network presented in the last paragraph according to experimental evidence reported in the literature. Subsequently, the simplified network is represented with the scheme proposed [21].

Finally, the models proposed by Crezee [22] and Zhang et al. [8] for the aldose hydrogenations are included to analyze the influence of the operating conditions on the yields of ethylene glycol, hexitols, and methane.

Simplified reaction network

The catalytic conversion of glucose to ethylene glycol can be represented according to the simplified scheme proposed [21]. This mathematical approach is based on a reaction network that only considers seven chemical reactions. For this purpose, the following simplifications are proposed by the cited authors for the global network:

A. Chemical reactions with yields below 7% are neglected. B. The production of mannitol and sorbitol are not distinguished because both reactions have a similar overall effect in the chemical process. C. Thermal side reactions are considered as two overall reactions for glucose and glycolaldehyde. D. Erythrose is considerably more reactive for the retro-aldol condensation than glucose.

In accordance with these simplifications, the global reaction network can be represented according to the scheme shown in Figure 1 and Table 1: The glycolaldehyde accumulation is approximately zero due to its high hydrogenation rate. In addition, the erythrose concentration is also considered quite low due to the minimum erythritol production observed experimentally by Zhao et al. [21]. According to the hypothesis D discussed above, each erythrose molecule produced per mole of glucose will transform into two molecules of glycolaldehyde instantaneously. Hence, the reaction 4 can also be omitted by considering that 3 molecules of glycolaldehyde are produced after the retro-aldol condensation of glucose.

Kinetic model N°1: Power law models

The kinetic model proposed by Zhao et al. [21] defines the rate of the reaction 𝑖 according to a power law model based on an apparent kinetic constant. This numerical approach is defined for the reactions 1, 3, 6 and 7:

reactant reactant exp i a i i i E r A c k c RT       (1)
Each power law corresponds to a pseudo-first order for the organic compound except for 𝑟 6 . Indeed, a different definition is required for the glycolaldehyde side reactions due to the high sensitivity of the ethylene glycol yield to the glycolaldehyde concentration. This condition is implemented in the mathematical model by considering a pseudo-second reaction order for the glycolaldehyde side reactions (𝑟 6 ).

The pre-exponential factor of the retro-aldol condensation (𝐴 1 ) can be modified to take into account the influence of the ammonium metatungstate concentration (CAMT) on the reaction rate.

The nonlinear increase of this rate with the catalyst amount can be represented with a modification of the value of the apparent kinetic constant [23]. In this manner, the term 𝐴 1 is substituted by the product of a corrected pre-exponential factor (𝐴 1 * ) and the AMT concentration at the power of 0.257:

* 0.257 11 AMT A A C  (2)
The value of the corrected pre-exponential factor is estimated in this study by considering the experimental tests carried out by Zhao et al. [21]. The results obtained in a batch reactor of 20 mL with 0.06 grams of AMT (CAMT = 1.01 mol•m -3 ) allowed establishing the value of 𝐴 1 * as 5.02•10 17 min -1 •(mol•m -3 ) -0.257 . Hereafter, the apparent kinetic constants of the reactions 1, 3, 6 and 7 are estimated according to the following parameters: Table 2. Parameters of the power law models [21], [23].

Kinetic model N°2: Langmuir Hinshelwood Hougen Watson (LHHW) equations

Three phases are involved in the process: Hydrogen, aqueous phase and solid catalyst. Hence, the hydrogenation reactions may be limited by mass-transfer. Glucose and glycolaldehyde react with adsorbed hydrogen at the surface of the supported metal. The hydrogen transport towards the liquid phase must be considered. The hydrogen solubility in aqueous solutions is quite low since D-glucose solutions are relatively viscous [22]. These issues affect the mass-transfer rate but are overcome by setting an impeller speed that is high enough to disperse a sufficient amount of gas in the liquid phase and increase the interfacial area with a greater bubble breakage. In this manner, the reaction rate becomes independent of the agitation regime and the first resistance is suppressed [8], [24]. Similarly, the liquid-solid mass transport effect shows efficiency above 95%, which indicates the possibility of neglecting the limitations on the transport towards the surface of the solid dispersed in the aqueous phase [22].

Secondly, the internal diffusion effects establish an important limitation only when the glucose concentrations are above 50% wt. or when the selected temperature implies a regime transition [22]. On the one hand, high glucose concentrations lead to a partial ionization of the molecule and therefore a low interaction with the solid catalyst. On the other hand, the temperature of the solution determines if the diffusive transport constitutes the controlling regime of the chemical reaction instead of the kinetics [START_REF] Červený | Catalytic hydrogenation[END_REF]. The chemical control condition was evidenced in previous studies that established that the activation energies of the hydrogenation reactions are greater than the activation energies of diffusion in liquids (12-21 kJ•mol -1 ) [7], [21]. However, the slopes of the Arrhenius's plots change due to a shift from a chemical control in the low-temperature region to a diffusion control in the high-temperature region [22]. This transition reduces the activation energies of the heterogeneous reactions. For this reason, it is necessary to consider different kinetic parameters for the hydrogenation reactions.

Hydrogenation can be described with pseudo-first order equations for glucose if the chemical process is carried out in the chemical control regime with an important hydrogen excess. For instance, Verma and Gehlawat [26] and Zhao et al. [21] described the kinetics of the D-glucose hydrogenation with nickel and ruthenium according to first and second order rate-laws based on apparent constants. However, this formulation cannot be envisaged for a wide temperature span due to the regime shift of the heterogeneous reactions.

Kinetic model formulation

A new definition was proposed for the hydrogenation rate-laws according to a Langmuir-Hinshelwood-Hougen-Watson (LHHW) equation [22]. This formulation provides an explicit way to describe the interaction of the reacting components with the catalytic surface [START_REF] Froment | Chemical reactor analysis and design[END_REF].

For this purpose, it is compulsory to take into account the predominance of each phenomenon associated with the mass transport and kinetics in the chemical reaction. In accordance with this statement, three kinetic models can be evaluated for the D-glucose hydrogenation [22]:

A. Non-competitive adsorption of hydrogen and D-glucose at different catalyst sites. The gas is adsorbed in a dissociative or molecular condition. B. Competitive adsorption of hydrogen and D-glucose at different catalyst sites. C. Competitive adsorption of dissociatively chemisorbed hydrogen and D-glucose.

The hydrogenation reactions are considered to be rate determining. Therefore, the corresponding overall rate equation of each kinetic model must also depend on the hydrogen pressure (𝑃 𝐻 ) [START_REF] Froment | Chemical reactor analysis and design[END_REF]. In accordance with this statement, the glucose hydrogenation can be described according to the following expression:

  22 1 H Glu Glu H n Glu Glu K K C P rk KC   (1)
The models are based on a linear dependency of the kinetic constant (𝑘 2 ) along with the adsorption constant of the glucose (𝐾 𝐺𝑙𝑢 ) and hydrogen (𝐾 𝐻 ). These alternatives are based on similar assumptions since they consider a negligible adsorption of solvent and products and do not make a distinction for the hexitols production. However, they differ from each other due to the exponent 𝑛 of their adsorption terms in the LHHW equation. The value of this variable is one, two or three for the alternatives A to C respectively.

The statistical analysis of these kinetic models established that all of them could be considered for the description of the hydrogenation reactions. Nevertheless, the alternative B seems unlikely since it considers a molecular hydrogen adsorption instead of a dissociative chemisorption [22].

In accordance with this statement, the exponent of the adsorption term can be one or three for the glucose hydrogenation.

The LHHW kinetic model shown in Equation 3 has been developed with an alternative formulation of the atomic hydrogen noncompetitive adsorption with the aldose adsorption [8]. By following this assumption, the reaction rates can be expressed with the following adaptations of the model A:

  ' 22 2 1 1 Glu Glu rut HH Glu Glu HH K C W KP rk K C V KP       (2)   ' 55 2 1 1 GA GA rut HH GA GA HH K C W KP rk K C V KP       (5)
[8], [START_REF] Zhang | [END_REF] established that the external diffusion limitations, associated with the gas-liquid and liquid-solid mass transfer, can be safely eliminated when the reaction was operated at a stirring speed above 900 rpm, catalyst loading below 0.6 g and catalyst particle size smaller than 74 mm (>200 mesh). In accordance with this statement, the experimental data are obtained under operating conditions that allow considering a low influence of the external mass transfer limitations. Moreover, the experimental analyses carried out by Zhang et al. [8] established that the hydrogen adsorption term in the denominator is negligible compared to the glucose adsorption term due to its considerably low concentration in the aqueous solution. This fact allows neglecting the factor associated with the hydrogen adsorption (1 + √𝐾 𝐻 𝑃 𝐻 ) and combining 𝐾 𝐻 with 𝑘 2 ′ or 𝑘 5 ′ .

Effects of the ammonium metatungstate on the heterogeneous catalyst

Some complementary modifications of the reaction rate-laws are proposed to take into account the suppressing effects exerted by the homogeneous catalyst on the aldose hydrogenations [8]. Firstly, it is necessary to consider the AMT adsorption on the surface of the ruthenium, which constitutes a catalyst poisoning effect. This can be achieved by including a term in the denominators of the rate-laws. For this purpose, Equations 4 and 5 must also have an additional adsorption factor based on the AMT concentration (𝐶 𝐴𝑀𝑇 ) and its adsorption equilibrium constant (𝐾 𝐴𝑀𝑇 ). Secondly, the binding characteristics of the ammonium metatungstate represent the formation of different complexes due to the interaction with the dissolved aldoses. For this reason, some additional assumptions can be considered to modify the LHHW rate-laws [8]:

 The complex formed by coordination of AMT with glucose or glycolaldehyde in solution can also be adsorbed on the Ru/AC surface.  Free or adsorbed AMT binds glucose molecules gradually until a maximum of four. On the contrary, glycolaldehyde is adsorbed with a stoichiometric ratio of 1/1.  The binding with AMT does not affect the surface reaction rate between the adsorbed aldose and hydrogen.

These assumptions define another term based on the adsorption terms of the rate-laws. These factors are determined by the equilibrium constants of the complexes formed with glucose (𝐾 𝐺-𝐴𝑀𝑇 ) and glycolaldehyde (𝐾 𝐺𝐴-𝐴𝑀𝑇 ). Henceforth, Equations 4 and 5 can be rewritten in the following manner:

' 3 10 Glu Glu C C  (6) '3 10 GA GA CC  (7) '3 10 AMT AMT CC  (8)   ' ' 4 ' 22 ' ' 4 ' ' () 1
( )

H H Glu AMT AMT Glu Glu Glu rut Glu AMT AMT Glu Glu Glu AMT AMT K P K C C K C W rk K C C K C K C V               (9)   ' ' ' 55 ' ' ' ' 1 H H GA AMT AMT GA GA GA rut GA AMT AMT GA GA GA AMT AMT K P K C C K C W rk K C C K C K C V               (10)
Each adsorption equilibrium and kinetic constant can be calculated according to a power law whose parameters are listed in Table 3: Table 3. Parameters for the equilibrium adsorption constants in the presence of AMT [8].

Constant

Compound Pre-exponential factor Activation or adsorption energy The association of this kinetic model with the power law equations requires the recalculation of the data listed in Table 3 according to the temperature range of the chemical process (180°C-240°C). For this purpose, the operation of a semi-batch reactor is represented with a mathematical model in order to estimate the values of the kinetic parameters of Equations 9 and 10 with the experimental data reported [5], [21].

(J•mol -1 ) G K Glucose 1.56•10 -3 -16500 GA K Glycolaldehyde 2.

Mathematical model of a semi-continuous reactor

The adaptation of the LHHW rate-laws in the global kinetic model is achieved with a numerical scheme based on an error minimization. This procedure is defined according to the root-meansquare error associated with the ethylene glycol and hexitols yields estimated with the constants of the LHHW model and the experimental data. The sum of the individual errors is minimized through an iterating scheme that adjusted the values of the constants listed in Table 3.

Firstly, this section describes the ordinary differential equations system that is envisaged to represent the experimental tests carried out in a semi-continuous reactor. The initial conditions of this system are established according to the experimental tests performed [21]. Then, the error minimization method is also explained in order to describe the scheme considered in this study for the parameters estimation.

Semi-continuous reactor

Glucose is mainly converted into hexitols (sorbitol and mannitol) when the chemical reactions are carried out in a batch reactor. Previous studies established that the ethylene glycol and hexitols yields are up to 7.0% and 39.5% respectively [5]. On the contrary, the EG yield reached an overall yield of 50.2% in a semi-continuous operation. This result indicates the necessity of a permanent feed at low flowrates that maintains low glucose concentrations to deplete the side reactions (2, 3 and 6). In addition, glucose readily undergoes a hydrogenation reaction to form stable sorbitol at temperatures as low as 80°C in the presence of Ru/AC catalyst [5].

Besides the selectivity increase generated by the glucose feed at high temperatures, a semibatch reactor also leads to higher yields due to the reactant concentration control. These operating conditions of the semi-batch reactor also lead to a reduction of the mass transfer limitations by reducing the viscosity of the liquid phase and generating glucose concentrations below 50 % wt. In the same manner, mass-transfer efficiency is enhanced with a high agitation speed (1000 rpm) and hydrogen pressures over 40 bar. In accordance with these statements, the numerical scheme of this study considers the reactor that is shown in Figure 2.

The system is modeled according to the perfectly stirred reactor approach proposed by Crezee [22], Zhang et al. [8] and Zhao et al. [21]. Therefore, the main aspects associated with the velocity profiles of the internal flow and the segregation levels of the heterogeneous catalyst are not considered in the computational model.

Ordinary differential equation system

The linear approach proposed in this study was modified in this study in order to describe the operation of a semi-continuous reactor according to an ordinary differential equation system. For this purpose, we defined a carbon mass balance for each compound by considering the rate-laws proposed for the 7 reactions of the simplified scheme. Firstly, the molar balance of glucose in a perfectly stirred semi-continuous reactor is defined in Equations 11 and 12. These expressions consider a constant feed rate of this saccharide (𝑁 𝐺𝑙𝑢 𝐹𝑒𝑒𝑑 ) and the rates of reactions 1 to 3.

3 Feed Glu Glu 1 i i dN N rV dt      (11) Feed Glu Glu0 N FC  (12)

5

Results and discussion

Kinetic parameters of the hydrogenation rate-laws

The kinetic model based on apparent constants that was proposed by Zhao et al. [21] determined that the regime transition of the controlling phenomenon occurs at different temperatures for both aldose hydrogenations. The glycolaldehyde hydrogenation attains the diffusion controlling regime when the temperature is above 100 °C whereas the glucose reaction requires a temperature over 160 °C to reach this level. This fact demands the adaptation of the rate-laws proposed by Zhang et al. [23] for a higher temperature range (180-240°C).

Parameters calculation

The kinetic parameters listed in Table 3 are recalculated with a numerical model based on the minimization of the least squares error. For this purpose, an objective function defines a global error as the sum of the yield errors of the main products: ethylene glycol, hexitols, and methane. The experimental data determined by Zhang et al. [5], [21] for a set of tests in a semi-batch reactor are compared with the yields predicted with the ordinary differential equations system. In this manner, the algorithm corrected the values of the kinetic parameters of Equations 9 and 10 by minimizing the sum of the squares of the errors. The values reported by Zhang et al. [8] are considered for the first iteration except for the activation energy and the pre-exponential factor of the product 𝑘 2 ′ 𝐾 𝐻 . These parameters are initially estimated according to the apparent constant proposed by Zhao et al. [21] for the constant 𝑘 2 in Equation 1. In this manner, the calculation considered that the true activation energy of the glucose hydrogenation is equal to twice its apparent activation energy [START_REF] Fogler | Elements of Chemical Reaction Engineering[END_REF]. Thereafter, the objective function is minimized by using the fminsearchbnd routine of MATLAB for nonlinear functions. This code allows considering bound constraints in the optimization engine used by fminsearch. For this purpose, a variable transformation is implemented through a wrapper function around the objective function. In this manner, each kinetic parameter of the LHHW equations is determined according to the thermodynamical aspects discussed by Zhang et al. [8]. The obtained values are reported in Table 4 The comparison of the kinetic data shown in Table 3 and Table 4 indicates a considerable diminution in the estimation of the glucose hydrogenation parameters leading to greater values for the parameter 𝑘 2 ′ 𝐾 𝐻 than those of the glucose adsorption constants. Moreover, the activation energy value of the kinetic factor of the LHHW equation approaches the activation energy value of diffusion in liquids. Thus, it is possible to establish that this chemical reaction is submitted to a diffusion controlling regime at the operating conditions of the one-pot reactor. On the contrary, the values of the glycolaldehyde hydrogenation do not indicate an important variation since the new estimation does not envisage a regime transition.

The kinetic model obtained by the modification of the LHHW equations is used in this section to analyze the influence of some key operating parameters on the selectivity of the chemical process. For this purpose, a sensitivity analysis is proposed in this study for the yields of ethylene glycol, hexitols, and methane. These results are determined by the following equations:

  0 3 EG EG Glu Glu VC Y FC VC   (13) 0 Hex Hex Glu Glu VC Y FC VC   (14)
 

0 6 Gas Gas Glu Glu VC Y FC VC   (15)
The parameter fitting is evaluated through a comparison of the yields predicted with the adjusted kinetic model and experimental data. This analysis is performed with the values listed in Table 5 for a set of semi-batch runs. Each experimental test was carried out for 40 minutes at various temperatures with 0.3 g of 4% Ru/AC. Then, the fitting parameters established a mean error equal to 5.8% on the predicted yields with regard to the experimental results.

Comparison with the power-law model

The behavior of the hydrogenation reactions that is predicted with the adjusted LHHW equations is compared with the experimental data obtained by Zhang et al. [5] for temperatures under the diffusion regime. Firstly, the results obtained in a batch reactor preloaded with glucose, 0.06 g of AMT and 0.3 g of 4% Ru/AC are compared with the concentration profile predicted by the adjusted LHHW equation (Figure 3A). The comparative analysis validates the adjustment of the kinetic model after establishing a root-mean squared error equal to 4•10 -4 mol•L -1 (2.3%) for an initial glucose concentration of 0.0175 mol•L -1 . Additionally, the concentration profiles predicted for glycolaldehyde hydrogenation have a good agreement with the experimental data obtained in the same batch reactor (Figure 3B). The comparison of the numerical model with the results obtained in a test carried out at 100 °C with an initial glycolaldehyde concentration of 8.5•10 -3 mol•L -1 indicates a root-mean squared error equal to 6.9•10 -4 mol•L -1 (8.1%). Similarly, a test performed at 120 °C with an initial concentration of 6•10 -3 mol•L -1 defines an error equal to 4.3•10 -4 mol•L -1 (7.2%) for the same analysis.

Influence of the reaction parameters

The variables that are considered in this section have a direct effect on the aldose hydrogenation rates. Nonetheless, only the glycolaldehyde hydrogenation belongs to the desired reaction pathway. Therefore, it is necessary to establish the appropriate level of these variables in order to enhance this reaction without increasing the undesired glucose hydrogenation considerably.

Previously, the influence of the reaction temperature and the glucose solution feed rate on the ethylene glycol yield has been analyzed [21]. Thus, we propose to study the following variables:

 Hydrogen partial pressure  Mass of homogeneous catalyst (ammonium metatungstate)  Mass of heterogeneous catalyst (ruthenium/activated carbon)

The simulation results show the final yields obtained in a semibatch reactor. The comparative study considers the results obtained after 120 minutes at three different temperatures: 180 °C, 220 °C and 240 °C. Additionally, the solution feed rate and the initial reactor volume are set to values of 0.667 mL•min -1 and 20 mL respectively. The spans of these sensitivity analyses are defined to describe the behavior of the reaction yields without generating significant variations on the dispersion and turbulence phenomena due to an excessive catalyst addition for the reaction volume. Thereafter, the simulation results of this section are obtained by setting two of the variables of this sensitivity analysis to a fixed value for the variation of the third one.

Hydrogen partial pressure

The hydrogen concentration in the aqueous solution is determined by its solubility. Thus, the influence of the hydrogen pressure can be described by Henry's law. The increase of the partial pressure of this gas defines a proportional augmentation of the hydrogenation rates [START_REF] Mishra | [END_REF]. Moreover, previous analyses varied the pressure from 40 to 75 bar and evidenced a linear dependence of this reaction rate on the pressure [22]. This fact allows establishing that this variable has a first-order behavior with respect to hydrogen. However, sugars hydrogenation over ruthenium is often described by zero order kinetics at high hydrogen pressures [30].

The hydrogenation rates determine the concentrations in the chemical reactor during the injection and reaction time. Therefore, the variations of the partial hydrogen pressure are directly associated with the selectivity of the semi-continuous process. Figure 4 describes the evolution of the concentrations predicted by the kinetic model for reactions at 35 and 60 bar. Conditions: AMT (0.06 g), Ru/AC 4% (0.3 g), 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 120 minutes.

Figure 4A and 4B indicate an important influence of the hydrogen pressure on the glucose and glycolaldehyde concentrations since these variables decrease when the reactor is pressurized at higher levels. This diminution is significant for both compounds when the reaction is carried out at 180 °C due to a low glucose conversion during the initial stages of the chemical process. Nevertheless, higher temperatures only define a significant reduction for the aldehyde. This result allows establishing how the ethylene glycol yield is reduced for excessive pressurization levels in the semibatch reactor. In the same manner, the behaviors observed for glucose and glycolaldehyde differ considerably from each other. On the one hand, Figure 4A and 4B indicate small variations of the glucose concentration during the reaction time. On the other hand, the glycolaldehyde concentrations show a significant time dependence since they reach peak values that are followed by a decreasing tendency at short reaction times. Evidently, low temperatures require longer times to reach the maximum glycolaldehyde concentrations due to a reduced cleavage rate in the glucose C-C bonds.

Moreover, the depletion of the glycolaldehyde concentration determines the transient behavior of the ethylene glycol concentration. The synthesis and hydrogenolysis of this diol are affected by high pressures when the reaction is carried out at 260 °C. During the first 60 minutes, the ethylene glycol production is increased by a higher presence of solubilized hydrogen in the liquid phase.

Thereafter, the concentration of this compound diminishes continuously due to the gas production within the reactor. In consequence, a limit reaction time must be considered for high reaction temperatures. This behavior is also expected at 180 °C and 220 °C but longer times are required to reach the maximum ethylene glycol concentration. These results allow establishing the most appropriate residence times when other operating conditions of the chemical reactor have already been determined.

Furthermore, the hydrogen pressure must be at least 20 bar to reduce the formation of byproducts such as humins or polymers [5]. However, it is recommended to pressurize the gas until 50 or 60 bar in order to obtain high ethylene glycol yields and a low gas production during the synthesis process. For this case study, Figure 5 describes the yields obtained for ethylene glycol, hexitols and methane at the temperatures that are considered above. Firstly, the results of the kinetic model at 260 °C (Figure 5A) indicate show an increase of the methane production when the hydrogen pressure is augmented. For pressures over 54 bar the gas production can achieve yields over 9% due to a high ethylene glycol decomposition. Therefore, the reactor operation above this pressure level is not recommended since it will be restrained by short injection periods. On the contrary, Figure 5C shows that the influence of the hydrogen pressure is negligible at 180 °C due to the high selectivity to hexitols. Finally, an intermediate behavior is observed at 220 °C in Figure 5B, which shows the increase in the production of the three compounds as the pressure increases. For this particular case, the operating pressure is defined by the augmentation of the three yields. In fact, the increases of the ethylene glycol yields stabilize when the pressure is over 47.5 bar whereas the increase of the by-products yields remains. This fact implies that pressures between 41 and 48 bar should be considered in order to reduce the hexitols and methane production and avoid an unnecessary increase of the operating costs. Conditions: AMT (0.06 g), Ru/AC 4% (0.3 g), 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 60 minutes.

Mass of homogeneous catalyst (ammonium metatungstate)

The influence of the homogeneous catalyst concentration on the ethylene glycol yield of this chemical process is also considered in this study through the variation of the amount of this catalyst that is preloaded in the reactor. Figure 6 shows the evolution of the concentrations of the organic compounds in the liquid phase when two different amounts of ammonium metatungstate are charged. Firstly, Figure 6A and 6B show that the glucose concentration decreases considerably due to the preload of higher concentrations of this catalyst if the reaction temperature is below 200 °C. On the contrary, the glycolaldehyde concentration shows an increase at all the reaction temperatures when the preloaded mass of AMT is 0.10 g instead of 0.02 g. The enhancement of the ethylene glycol production that is associated with the addition of greater amounts of AMT is inferior for reactions at higher temperatures. For instance, Figure 7 shows that reactions carried out at 220 °C have an increase of the glycolaldehyde concentration with higher AMT concentrations. Nonetheless, this augmentation does not represent a significant increase in the diol concentration when the equipment is operated at 260 °C. This result is due to the thermal side reactions of this compound, which are also enhanced by the heating of the chemical reactor. This aspect is more evident when the equipment is operated at this temperature since the aldehyde concentration is considerably lower than that of 220 °C. Hence, the glycol concentration is not significantly increased either.

Furthermore, Figures 6C and6D show that the addition of AMT extends the limit reaction time at which the diol decomposition occurs, from 60 to 80 minutes, when the reaction is carried out at 260 °C. For this reason, this parameter must be mainly associated with the reaction temperature and the influence of the concentration of this catalyst on the reaction yields must be established according to the operating temperature (Figure 7). For example, Figure 7A shows that the hexitols production is reduced from 2.2% to only 0.6% with the addition of 0.06 g of AMT instead of 0.02 g at 260 °C. On the contrary, an operation at 220 °C reduces the hexitols yield from 21.6% to 7.9% with the same catalyst addition (Figure 7B). Finally, the evaluation of the yields obtained at 180 °C (Figure 7C) indicates that the increase of the ethylene glycol concentration is significant but insufficient since the change from 0.02 g to 0.10 g of AMT only represents an augmentation of 8.6% in the diol yield. These results allow concluding that the AMT concentration can be an important operating parameter only at intermediate temperatures (200-240 °C). Conditions: Ru/AC 4% (0.3 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 60 minutes.

Mass of heterogeneous catalyst (ruthenium/activated carbon)

The influence of the heterogeneous catalyst is also considered in this study through the variation of the amount of Ru/AC initially preloaded in the chemical reactor. The effects of a variation of this operating parameter are correlated to the glucose and glycolaldehyde hydrogenation rates. This influence can be analyzed with the kinetic model if the variation of the catalyst concentration does not constitute significant changes in the solid's dispersibility or the turbulence regime of the reactive mixture. In accordance with this statement, Figure 8 illustrates the concentrations of aldoses, hexitols and ethylene glycol that are obtained by varying the mass of Ru/AC between 0.1 and 0.5 g. The addition of higher concentrations of the heterogeneous catalyst has the opposite behavior to that observed with the increase of the AMT concentration (Figure 6Erreur ! Source du renvoi introuvable.). A decrease of the glycolaldehyde and aldehyde concentrations is observed in Figure 8A and 8B when the initial amount of Ru/AC is increased. Moreover, the simulation results indicate how the Ru/AC concentration defines the influence of the reaction temperature. On the one hand, low Ru/AC concentrations determine the glycolaldehyde and glucose concentrations according to the reaction temperature. On the other hand, high Ru/AC concentrations do not generate significant changes in the concentrations of both aldoses when the reactor is heated in spite of the enhancement of the retro-aldol condensation. This result is due to the immediate hydrogenation of glucose and glycolaldehyde that overcomes the effect associated with the temperature change. This fact illustrates how the mass of Ru/AC preloaded in the equipment determines the transient concentrations of glycolaldehyde and glucose.

Figure 8A shows high accumulations of glycolaldehyde and glucose in the aqueous solution when only 0.1 grams of Ru/AC are preloaded in the system. Figure 8C indicates that the maximum ethylene glycol concentration that is obtained in the process corresponds to only 0.42 mol•L -1 (220 °C). This result allows establishing that a minimum amount of Ru/AC must be charged in the equipment in order to avoid an excessive glycolaldehyde accumulation. Indeed, an operation without the minimum amount of Ru/AC leads to greater yields of the thermal side reactions.

Furthermore, Figure 8Erreur ! Source du renvoi introuvable.D illustrates the main issue that can be associated with an excess of the heterogeneous catalyst in the mixture. Despite the small variations of the glycolaldehyde and glucose concentrations during the reaction with an excess of Ru/AC, the ethylene glycol production has a transient behavior that depends on the Ru/AC concentration. Indeed, an operation at 260 °C enhances the diol production but increases its degradation rate after the limit reaction time. In addition, these operating conditions also define the greatest decrease of the ethylene glycol concentration after 70 minutes. On the contrary, an operation at 180 °C has a low time dependence due to the high selectivity to hexitols. Finally, the operation at 220 °C shows that the maximum ethylene glycol concentration that is evidenced in this study (0.93 mol•L -1 in Figure 8D) is obtained with the addition of 0.5 g of Ru/AC.

The behaviors observed with the addition of an insufficient and excessive amount of Ru/AC allow concluding that it is possible to define a set of operating conditions that maximize the ethylene glycol yield without compromising the reaction time. Figure 9 shows the reaction yields that can be obtained at different temperatures and Ru/AC concentrations. This chart shows that all the reaction yields increase by the addition of more Ru/AC to the mixture. For this reason, the mass of heterogeneous catalyst that is charged into the chemical reactor must be determined according to a previous setting of the reaction time and temperature. Firstly, Figure 9A shows that an operation at 260 °C that lasts 60 minutes must not be carried out with more than 0.1 g of Ru/AC because the methane yield can be greater than 5%. Secondly, Figure 9B shows that an operation at 220 °C allows increasing the mass of this catalyst until 0.4 g before reaching a hexitols yield over 10%. Finally, Figure 9C shows that an operation at 180 °C must not consider high Ru/AC concentrations because this variation only enhances the hexitols production. Conditions: AMT (0.06 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 60 minutes.

Conclusions

The development of a combined kinetic model for the conversion of glucose into ethylene glycol allows analyzing the influence of its operating conditions on the reaction yields. For this purpose, the kinetic parameters of the hydrogenations rate-laws were adjusted to the temperature range that enhances the ethylene glycol production (180 °C -240 °C). Subsequently, this numerical model was used to evaluate the glucose conversion and selectivity in a semibatch reactor. This analysis identified the relation between the operating conditions and the transient concentration profiles of the reagent and the main products of this chemical process.

Firstly, the kinetic model allows establishing that the reaction temperature defines not only the selectivity but also the time in which the maximum ethylene glycol concentration is obtained in the equipment. In accordance with these results, the operating temperature that is recommended for high selectivity to ethylene glycol must be set at 220 °C approximately. The operation at this temperature defines a period during which a maximum ethylene glycol yield can be obtained along with a low hexitols and methane production. This enhancement can be improved by adding greater concentrations of the homogeneous catalyst in the mixture. However, the selectivity of the process can be negatively affected by operations carried out for periods longer than this critical time.

Secondly, other parameters such as the hydrogen partial pressure and the Ru/AC concentration have an impact on the maximum ethylene glycol concentration that can be reached in the semibatch reactor and the maximum reaction time. For instance, an increase in the reaction rates can be achieved with higher hydrogen and Ru/AC concentrations. Nonetheless, glucose hydrogenation is enhanced in a greater proportion than that of glycolaldehyde hydrogenation. Therefore, a maximum level must be defined for these variables in order to avoid an excessive conversion of glucose into hexitols. These limit concentrations are determined by the other operating conditions of the chemical reactor.

In conclusion, the combination of the two distinct approaches in the kinetic model is an important step forward in the development of a more accurate model. This integration allows analyzing the influence of the hydrogen pressure and the catalyst concentrations on the reaction selectivity. Furthermore, a detailed description of the flow characteristics during the agitation process can contribute to the definition of the residence time through the analysis of phenomena such as liquid re-circulation or parabolic gas holdup. For this purpose, further stages of this study will focus on the flow characterization and the influence of the operating conditions on the segregation levels of the heterogeneous catalyst. In this manner, the model defines the main aspects that must be considered during the design of a unit that benefits from the conversion of non-alimentary glucose for the production of ethylene glycol. Table 2. Parameters of the power law models [21], [23].

Symbols used

i A [mol -1 • L -1 •s -1 ] or [s -1 ]
Table 3. Parameters for the equilibrium adsorption constants in the presence of AMT [8].

Table 4. Estimated parameters for the equilibrium adsorption constants in presence of AMT (temperatures between 180 °C and 260 °C).

Table 5. Comparison of the yields estimated with the adapted kinetic model and the experimental data [21].

Entry for the Table of Contents

Figure 1 .

 1 Figure 1. Simplified reaction network [21].

Figure 2 .

 2 Figure 2. Scheme of the semi-batch reactor and the solution injection system.

Figure 3 .

 3 Figure 3. Comparison of the aldose concentrations predicted with the hydrogenation models in a batch reactor.

Figure 4 .

 4 Figure 4. Influence of the hydrogen pressure on the concentrations in the liquid phase.

Figure 5 .

 5 Figure 5. Influence of the hydrogen pressure on the reaction yields.

Figure 6 .

 6 Figure 6. Influence of the mass of AMT on the concentrations in the liquid phase. Conditions: Ru/AC 4% (0.3 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 120 minutes.

Figure 7 .

 7 Figure 7. Influence of the mass of AMT on the reaction yields.

Figure 8 .

 8 Figure 8. Influence of the mass of Ru/AC 4% on the concentrations in the liquid phase. Conditions: AMT (0.06 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 120 minutes.

Figure 9 .

 9 Figure 9. Influence of the mass of Ru/AC 4% on the reaction yields.

1 ARj

 1 Pre-exponential factor of the chemical reaction in the power law model * [min -1 •(mol•m -3 ) -0.257 ] Corrected pre-exponential factor of the retro-L -1 ]Concentration of the j compound in the chemical reactor' j C [mol•m -3 ]Concentration of the j compound in the chemical reactor 1•L•s -1 ] or [s -1 ] Kinetic constant of the reaction 𝑖 (Power-law equation) AMT K [m 3 •mol -1 ] AMT adsorption equilibrium constant on Ru/AC GA K [m 3 •mol -1 ] Glycolaldehyde adsorption equilibrium constant on Ru/AC GA AMT K  [m 6 •mol -2 ] GA-complexes adsorption equilibrium constant on Ru/AC Glu K [m 3 •mol -1 ] Glucose adsorption equilibrium constant on Ru/AC Glu AMT K  [m 15 •mol -5 ] Glucose-complexes adsorption equilibrium constant on Ru/AC [kJ•mol -1 •K] Ideal gas constant i r [mol•L -1 •s -1 ] Rate of the reaction 𝑖 T of by-products produced per mole of glucose y [-] Moles of by-products produced per mole of glycolaldehyde j Y [-] Yield of the compound j on the chemical reactions network Greek letters  [-] Reaction order of the chemical reaction in the power law equation , ji  [-] Stoichiometric coefficient of the compound 𝑗 in the reaction 𝑖 Compound (𝑗 = 𝐺𝑙𝑢, 𝐻𝑒𝑥, 𝐸𝑟𝑦, 𝐺𝐴, 𝐸𝐺, 𝐺-𝐵𝑦𝑃, 𝐺𝐴-𝐵𝑦𝑃, 𝐺𝑎𝑠) (Glucose, hexitols, erythrose, glycolaldehyde, ethylene glycol, glucose by-products, glycolaldehyde by-products, methane) [30] V. Sifontes, O. Oladele, K. Kordás, K. Eränen, J. Mikkola, D. Murzin, T. Salmi, J. Chem. Technol. Biotechnol. 2011, 86 (5), 658-668. DOI: 10.1002/jctb.2565.

Figure captions Figure 1 .

 captions1 Figure captionsFigure1. Simplified reaction network[21].

Figure 2 .

 2 Figure 2. Scheme of the semi-batch reactor and the solution injection system.

Figure 3 .

 3 Figure 3. Comparison of the aldose concentrations predicted with the hydrogenation models in a batch reactor.

Figure 4 .

 4 Figure 4. Influence of the hydrogen pressure on the concentrations in the liquid phase.Conditions: AMT (0.06 g), Ru/AC 4% (0.3 g), 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 120 minutes.

Figure 5 .

 5 Figure 5. Influence of the hydrogen pressure on the reaction yields.Conditions: AMT (0.06 g), Ru/AC 4% (0.3 g), 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 60 minutes.

Figure 6 .

 6 Figure 6. Influence of the mass of AMT on the concentrations in the liquid phase.Conditions: Ru/AC 4% (0.3 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 120 minutes.

Figure 7 .

 7 Figure 7. Influence of the mass of AMT on the reaction yields.Conditions: Ru/AC 4% (0.3 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 60 minutes.

Figure 8 .

 8 Figure 8. Influence of the mass of Ru/AC 4% on the concentrations in the liquid phase. Conditions: AMT (0.06 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 120 minutes.

Figure 9 .

 9 Figure 9. Influence of the mass of Ru/AC 4% on the reaction yields. Conditions: AMT (0.06 g), 50 bar H2, 10% glucose solution, initial volume: 20 mL, feed flowrate: 0.667 mL.min -1 , reaction time: 60 minutes.

Table 1 .

 1 T Simplified chemical reaction network of the catalytic conversion of glucose to ethylene glycol[21].

	Chemical reaction	Chemical equation
	1. Retro-aldol condensation of glucose to erythrose and glycolaldehyde	6 12 6 C H O	4 8 4 C H O C H O 2 4 2 
	2. Glucose hydrogenation to hexitols	C	H		O	H 	C	H	O
		6	12	6				2			6	14	6
	3. Glucose side reactions	C H O		xByP
		6 12 6
	4. Erythrose conversion to glycolaldehyde	C	H	O		2 C	H	O
		4	8		4					2	4	2
	5. Glycolaldehyde hydrogenation to ethylene glycol	C	H	O	H 	C	H	O
		2	4	2				2			2	6	2
	6. Glycolaldehyde side reactions	C	H	O		yByP
		2	4		2
	7. Methane production	C	H	O		3	H		2 CH		2	H	O
		2	6		2				2					4	2

Table 4 .

 4 : Estimated parameters for the equilibrium adsorption constants in presence of AMT (temperatures between 180 °C and 260 °C).

	Constant	Compound	Pre-exponential factor	Activation or adsorption energy (J•mol -1 )
	K	Glucose	5.91•10 -3	-18577
	G
	K	Glycolaldehyde	1.00•10 -12	-2699
	GA
	Glu AMT K 	Complexes formed with glucose	1.87•10 -11	-63167

Table 5 .

 5 Comparison of the yields estimated with the adapted kinetic model and the experimental data[21].

	Mass of AMT -Ru/AC (g)	Temperature (K)	H2 pressure (bar)	Glucose concentration in feed (% wt.)	Feed rate (mL•min -1 )	Experimental yields (%) EG Hex.	CH4	Kinetic model (%) EG Hex.	CH4	Absolute error (%) EG Hex.	CH4
	0.060 -	453	50	10	0.667	30.0	65.1	0.0	30.1	64.7	0.1	-0.1	0.4	-0.1
	0.3
	0.060 -	473	50	10	0.667	59.5	32.0	0.7	59.8	25.7	0.6	-0.3	6.3	0.1
	0.3

Table captions Table 1 .

 captions1 T Simplified chemical reaction network of the catalytic conversion of glucose to ethylene glycol[21].

Acknowledgement

Nicolas Godin, Isariebel Quesada and Laurence Llovera from the Industrial Company of Vegetable Matter (CIMV) and Idaline Chanteperdrix from the Laboratory of Coordination Chemistry (LCC) are gratefully acknowledged for their support and contribution to the development of this scientific study.

Type of Article:


Thanida Sananmuang

Denis Puthier

Catherine Nguyen

Kaj Chokeshaiusaha

Novel Classifier Orthologs Of Bovine And Human Oocytes Matured In Different Melatonin Environments

Keywords: Classifier orthologs, Melatonin, Oocytes, Cross-species RNA-Seq meta-analysis

It has been demonstrated that melatonin influences the developmental competence of both in vivo and in vitro matured oocytes. It modulates oocyte-specific gene expression patterns among mammalian species. Due to differences among study systems, the identification of the classifier orthologs-the homologous genes related among mammals that could universally categorize oocytes matured in environments with varied melatonin levels is still limitedly studied. To gain insight into such orthologs, cross-species transcription profiling meta-analysis of in vitro matured bovine oocytes and in vivo matured human oocytes in low and high melatonin environments was demonstrated in the current study. RNA-Seq data of bovine and human oocytes were retrieved from the Sequence Read Archive database and pre-processed. The used datasets of bovine oocytes obtained from culturing in the absence of melatonin and human oocytes from old patients were regarded as oocytes in the low melatonin environment (Low). Datasets from bovine oocytes cultured in 10-9 M melatonin and human oocytes from young patients were considered as oocytes in the high melatonin environment (High). Candidate orthologs differentially expressed between Low and High melatonin environments were selected by a linear model, and were further verified by Zero-inflated regression analysis. Support Vector Machine (SVM) was applied to determine the potentials of the verified orthologs as classifiers of melatonin environments. According to the acquired results, linear model analysis identified 284 candidate orthologs differentially expressed between Low and High melatonin environments. Among them, only 15 candidate orthologs were verified by Zero-inflated regression analysis (FDR ≤ 0.05). Utilization of the verified orthologs as classifiers in SVM resulted in the precise classification of oocyte learning datasets according to their melatonin environments (Misclassification rates < 0.18, area under curves > 0.9). In conclusion, the cross-species RNA-Seq meta-analysis to identify novel classifier orthologs of matured oocytes under different melatonin environments was successfully demonstrated in this study-delivering candidate orthologs for future studies at biological levels. Such verified orthologs might provide valuable evidence about melatonin sufficiency in target oocytes-by which, the decision on melatonin supplementation could be implied.

(ART) in both humans and animal species [START_REF] Wang | in vitro fertilization (IVF): A review of 3 decades of clinical innovation and technological advancement[END_REF]. in vitro oocyte maturation (IVM) is one of the multiple steps of IVP that the oocytes are recovered directly from the follicles and are cultured in a specific condition. Despite its advancement, incomparable in vitro and in vivo environment is still a major problem contributing to the success of IVM performance [START_REF] Zhang | Insufficiency of melatonin in follicular fluid is a reversible cause for advanced maternal age-related aneuploidy in oocytes[END_REF][START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF]. During in vivo oocyte development, several delicate substances and hormones are locally secreted in follicular fluid by ovarian and granulosa cells to nourish oocytes. This includes melatonin which was widely studied for its contribution to oocyte quality [START_REF] Zhang | Insufficiency of melatonin in follicular fluid is a reversible cause for advanced maternal age-related aneuploidy in oocytes[END_REF][START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF][START_REF] Do | Melatonin Supplementation During in vitro Maturation and Development Supports the Development of Porcine Embryos[END_REF][START_REF] Soto-Heras | Effects of melatonin on oocyte developmental competence and the role of melatonin receptor 1 in juvenile goats[END_REF]. Melatonin was most recognized for its scavenging function on reactive oxygen species (ROS). Since ROS play a major role in impairing oocyte development [START_REF] El-Raey | Evidence of melatonin synthesis in the cumulus oocyte complexes and its role in enhancing oocyte maturation in vitro in cattle[END_REF][START_REF] Jing | Melatonin levels in follicular fluid as markers for IVF outcomes and predicting ovarian reserve[END_REF], the lack of melatonin in conventional in vitro culture environments has contributed to ROSinduced poor oocyte quality in several mammals [START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF][START_REF] Lee | Melatonin influences the sonic hedgehog signaling pathway in porcine cumulus oocyte complexes[END_REF][START_REF] Igarashi | Oocyte aging underlies female reproductive aging: biological mechanisms and therapeutic strategies[END_REF].

Several factors also affect in vivo oocyte development. This includes reproductive aging-one of the most common causes of female infertility. It has been reported that reproductive aging induces a variety of physiological aberrations withering oocyte quality [START_REF] Igarashi | Oocyte aging underlies female reproductive aging: biological mechanisms and therapeutic strategies[END_REF][START_REF] Malhi | A bovine model to study reproductive aging[END_REF][START_REF] Malhi | Oocyte developmental competence in a bovine model of reproductive aging[END_REF][START_REF] Thouas | Effect of Female Age on Mouse Oocyte Developmental Competence Following Mitochondrial Injury1[END_REF]. One of them is melatonin deprivation, caused by reduced melatonin production from aging granulosa and ovarian cells. Melatonin deprivation was hereby blamed as a major contributor to aging oocyte's incompetency [START_REF] Zhang | Insufficiency of melatonin in follicular fluid is a reversible cause for advanced maternal age-related aneuploidy in oocytes[END_REF][START_REF] Igarashi | Oocyte aging underlies female reproductive aging: biological mechanisms and therapeutic strategies[END_REF]. Aging also results in other diverse hormonal imbalances and assorted physiological changes [START_REF] Igarashi | Oocyte aging underlies female reproductive aging: biological mechanisms and therapeutic strategies[END_REF][START_REF] George | Fertility and age[END_REF]. All these factors represent complex negative effects during in vivo oocyte aging. Therefore, it has been difficult to exclusively determine the melatonin deprivation effects on the oocyte during in vivo development in such circumstances. By means of this, it was difficult to exclusively determine the melatonin deprivation effect on in vivo oocyte in such circumstances.

Interestingly, similar aging aberrance also occurred with in vitro cultured oocytes in a time-dependent manner. in vitro-aging phenomenon also included ROS-induced quality deterioration due to the absence of melatonin in the conventional culture environment [START_REF] Zhang | Insufficiency of melatonin in follicular fluid is a reversible cause for advanced maternal age-related aneuploidy in oocytes[END_REF][START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF][START_REF] Do | Melatonin Supplementation During in vitro Maturation and Development Supports the Development of Porcine Embryos[END_REF][START_REF] Soto-Heras | Effects of melatonin on oocyte developmental competence and the role of melatonin receptor 1 in juvenile goats[END_REF][START_REF] Igarashi | Oocyte aging underlies female reproductive aging: biological mechanisms and therapeutic strategies[END_REF][START_REF] Malhi | Oocyte developmental competence in a bovine model of reproductive aging[END_REF]. In agreement with the protective effects of in vivo melatonin produced by ovarian follicles, in vitro melatonin supplement was evidenced to restore the oocyte developmental competency in several species [START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF][START_REF] Do | Melatonin Supplementation During in vitro Maturation and Development Supports the Development of Porcine Embryos[END_REF][START_REF] Soto-Heras | Effects of melatonin on oocyte developmental competence and the role of melatonin receptor 1 in juvenile goats[END_REF][START_REF] Igarashi | Oocyte aging underlies female reproductive aging: biological mechanisms and therapeutic strategies[END_REF]. Unlike in vivo-aging condition divided effects of in vitro melatonin deprivation should be conveniently observed due to the simplicity of the culture environment.

Through this, analysis of crossed effects between in vitro and in vivo melatonin deprivation due to aging should hereby contribute a novel insight into the universal melatonin effects independent from other factors.

With comparable reproductive physiology and embryogenesis, cattle are regarded as excellent study models for human IVP [START_REF] Malhi | A bovine model to study reproductive aging[END_REF][START_REF] Hansen | Realizing the promise of IVF in cattle -An overview[END_REF]. Employing this, novel transcriptomic technologies including RNA sequencing (RNA-Seq) have been continuously applied with oocyte samples obtained from both species for studying their revised universal transcriptome's regulation. However, the application of standard RNA-Seq or bulk RNA-Seq with mammalian oocytes was not practical due to low oocyte numbers for cDNA library preparation [START_REF] Brayboy | The double-edged sword of the mammalian oocyte -advantages, drawbacks and approaches for basic and clinical analysis at the single cell level[END_REF][START_REF] Chambers | Transcriptomics and single-cell RNA-sequencing[END_REF]. Moreover, oocytes were heterogeneous populations, and should thus not represented with average results. Since cell-to-cell variations are common among isolated oocytes, loss of real features among them in bulk analysis could bring about considerable bias in the analysis [START_REF] Brayboy | The double-edged sword of the mammalian oocyte -advantages, drawbacks and approaches for basic and clinical analysis at the single cell level[END_REF].

Recently, low-input and single-cell RNA sequencing (RNA-Seq) methods have been implemented with oocytes acquired from various mammals [START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF][START_REF] Brayboy | The double-edged sword of the mammalian oocyte -advantages, drawbacks and approaches for basic and clinical analysis at the single cell level[END_REF][START_REF] Liu | Single-cell transcriptome sequencing reveals that cell division cycle 5-like protein is essential for porcine oocyte maturation[END_REF]. Since such technologies allow robust transcript profiling analysis down to a single cell, they enable researchers to explore biological variations among oocytes with reduced bias [START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF][START_REF] Chambers | Transcriptomics and single-cell RNA-sequencing[END_REF][START_REF] Chitwood | 131 RNA-Seq TRANSCRIPTOME PROFILING OF INDIVIDUAL RHESUS MACAQUE OOCYTES AND PRE-IMPLANTATION EMBRYOS[END_REF]. Interestingly, an increase of such data among a variety of mammalian species even allowed cross-species meta-analysis of oocytes. Since cross-species transcription profiling meta-analysis of oocytes could imply universal orthologs regulated among the species in the environment of interest [START_REF] Chokeshaiusaha | Construction of pooled oocyte expression profiles of rhesus monkey and mouse for concurrent meta-analyses[END_REF], implement of such analysis with RNA-Seq data of oocytes obtained from different in vitro and in vivo melatonin environments should imply us insight into such orthologs, as well.

In several transcription profiling studies, genomic classification or subtyping of unknown samples can be archived by various machine learning algorithms. Support vector machine (SVM) is a robust learning model widely implemented on a genome-wide scale to prioritize orthologs for sample classification [START_REF] Slawski | CMA -A comprehensive Bioconductor package for supervised classification with high dimensional data[END_REF][START_REF] Fernandes De Mello | Introduction to Support Vector Machines[END_REF]. Novel classifier orthologs were identified by SVM in several topics including aging, cancers, infections, and metabolic disorders [START_REF] Huang | Applications of support vector machine (SVM) learning in cancer genomics[END_REF][START_REF] Wang | Machine learning-based method for obesity risk evaluation using single-nucleotide polymorphisms derived from next-generation sequencing[END_REF][START_REF] Kumar | SVMRFE based approach for prediction of most discriminatory gene target for type II diabetes[END_REF][START_REF] Arabfard | Genome-wide prediction and prioritization of human aging genes by data fusion: A machine learning approach[END_REF]. Regardless of its flexibility, SVM application in oocyte study was still limitedly demonstrated.

It was noteworthy that the integration of SVM with cross-species meta-analysis of different melatonin environments could thus help us determining novel classifier orthologs for evaluating the melatonin sufficiency of oocytes. In other words, such cognition would preliminarily assess the beneficial effects of melatonin supplement on target oocytes presented in environments of interest.

To gain insight into universal classifier orthologs for melatonin environments, cross-species transcription profiling meta-analysis of in vitro matured bovine oocytes and in vivo matured human oocytes in low and high melatonin environments was performed in the current study. The pooled cDNA library based on human-bovine orthologs was constructed with bias correction. With restricted verification procedures and SVM classification, the universal classifier orthologs for melatonin environments of oocytes were successfully implied for the first time in this study.

Materials and methods revised

Sample datasets

RNA-Seq data of bovine and human oocytes were retrieved from the sequence read archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra) (Table 1). Both bovine and human RNA-Seq data were presented in NCBI with accession numbers GSE122738 [START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF] and GSE125300 [START_REF] Zhang | Advanced maternal age alters expression of maternal effect genes that are essential for human oocyte quality[END_REF] (https://www.ncbi.nlm.nih.gov/geo/).

The sources of bovine [START_REF] An | Melatonin supplementation during in vitro maturation of oocyte enhances subsequent development of bovine cloned embryos[END_REF] and human oocytes [START_REF] Zhang | Advanced maternal age alters expression of maternal effect genes that are essential for human oocyte quality[END_REF] utilized for data preparation were summarized in Table 1.

In brief, bovine cumulus-oocyte complexes (COCs) were collected from slaughterhouses and cultured in the presence (In_vitro_MII+Mela) or absence (In_vitro_MII) of 10 -9 M melatonin for 24 hours, respectively.

Metaphase II (MII) oocytes were isolated from COCs and divided into 6 samples-each containing 7 oocytes of their groups. On the contrary, each human MII oocyte sample was retrieved directly from single follicles of young (In_vivo_young_MII), or old (In_vivo_old_MII) patients by ultrasound-guided technique. In this study, the term-melatonin environment was used to infer the amount of melatonin in the environment surrounding the target oocyte samples. Through this, the bovine oocytes matured in the absence of melatonin supplement, and human oocytes acquired from old patients were regarded as oocytes in the low melatonin environment. On the other hand, bovine oocytes matured in the presence of melatonin supplement, and human oocytes acquired from young patients were considered as oocytes in the high melatonin environment.

Data pre-processing

Data pre-processing was performed with the process similar to that described in our previous study [START_REF] Chokeshaiusaha | Factor Analysis for Bicluster Acquisition (FABIA)revealed vincristine-sensitive transcript pattern of canine transmissible venereal tumors[END_REF]. In brief, all sequences were trimmed to remove contaminated adapter sequences and unqualified sequences (length < 25 nucleotides and mean Phred score < 25). Duplicated sequences removal, sequence alignments (>85%), and GC base-balance bias correction were accomplished prior to gene-level counting. In this study, the read counts were normalized by the trimmed mean of M-values (TMM). bovine-human orthologs were assigned with the previously reported procedure [START_REF] Chokeshaiusaha | Construction of pooled oocyte expression profiles of rhesus monkey and mouse for concurrent meta-analyses[END_REF]. Data quality control was carried out [START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF] along with cell-cycle effect removal to minimize the bias from within-cell-type heterogeneity [START_REF] Barron | Identifying and removing the cell-cycle effect from single-cell RNA-Sequencing data[END_REF]. Count read numbers of all samples were > 14 x10 -6 , and none of them was prone to apoptosis (Fraction of mitochondrial reads < 1.5%).

The mutual nearest neighbors (MNN) in ortholog expression space was identified for cross-species batch correction [START_REF] Shafer | Cross-Species Analysis of Single-Cell Transcriptomic Data[END_REF] of pooled cDNA library created from bovine and human oocyte count data using 'scater' package [START_REF] Haghverdi | Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors[END_REF].

Candidate ortholog selection for oocyte sample categorization based on melatonin environments revised

Candidate orthologs for oocyte sample categorization by melatonin environment were selected using procedures as described in our previous study [START_REF] Chokeshaiusaha | Construction of pooled oocyte expression profiles of rhesus monkey and mouse for concurrent meta-analyses[END_REF]. In brief, candidate bovine-human orthologs differentially expressed between oocytes in low and high melatonin environments were selected from 1,000 learning datasets using the linear model and empirical Bayes methods. In this study, we consider differentially expressed orthologs with importance value ≥ 2.5. Only orthologs declared significant both by intra-species and cross-species analyses would be selected as candidate orthologs. With expression values of the candidate orthologs, Partitioning around medoids (PAM) clustering (2-4 clusters) was performed, and clustering validation was determined by Connectivity, Dunn, and Silhouette scores.

Verification of the candidate orthologs

Concerning about unrealistic mean-variance relationship produced by the batch-corrected data, the candidate orthologs were verified by Zero-inflated regression analysis with proper model matrix using raw count data to avoid the bias from the batch correction process with 'MAST' package [START_REF] Finak | MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data[END_REF]. The same criteria of candidate ortholog selection was applied-by which the verified candidate orthologs must be differentially expressed between low and high melatonin environments both by intra-species and cross-species analyses (FDR ≤ 0.05).

Oocyte sample classification by Support Vector Machine

To determine the potential of the verified candidate orthologs as classifiers for melatonin environments, Support Vector Machine (SVM) with linear kernel classifier was performed. Different training data were generated by 5-fold cross-validation, Bootstrap, and Monte-Carlo-cross-validation strategies-1,000 iterations for each strategy. All procedures could be archived by class prediction functions using the 'CMA' package [START_REF] Slawski | CMA -A comprehensive Bioconductor package for supervised classification with high dimensional data[END_REF].

Data visualization

Heatmap illustrating candidate orthologs was drawn using the "ComplexHeatmap" package [START_REF] Gu | Complex heatmaps reveal patterns and correlations in multidimensional genomic data[END_REF]

Oocyte samples could be clustered by their melatonin environments using the candidate orthologs as features.

Analysis of bias-corrected data by linear model revealed 284 candidate orthologs differentially expressed between oocytes in low and high melatonin environments (importance value ≥ 2.5). According to PAM clustering results based on expressions of these candidate orthologs (Table 2 andFig

All verified candidate orthologs were significantly expressed in oocytes in the high melatonin environment.

Zero-inflated regression analysis revealed 158 orthologs differentially expressed between oocytes in low and high melatonin environments (FDR≤0.05). However, only 15 orthologs were found intersected with the candidate orthologs acquired from the previous selections-AAMP, ADM, ATP6V0A2, BCL7B, CYB561D2, ENDOV, F2RL1, KIAA0586, LPCAT4, NR1H3, PLEKHM1, SLC47A1, SYT11, RFT1, and MRPS18C (Fig 3). While all oocytes in the high melatonin environment expressed these orthologs higher than those in the low environment, only ADM, CYB561D2, ENDOV, KIAA0586, NR1H3, PLEKHM1, SLC47A1, SYT11, RFT1, and MRPS18C rendered substantial fold-change values (≥ 1.5). According to our review of their biological functions, such orthologs were associated with a variety of cell metabolisms (Table 3).

Verified candidate orthologs were potential classifiers for the melatonin environments of oocytes.

The verified candidate orthologs were used as classifier orthologs in SVM performance. Consistent results were acquired among learning datasets-by which all of them rendered Misclassification rates < 0.18, and their ROCs had Area Under Curves (AUCs) > 0.9 (Fig 4).

Discussion

A cross-species transcription profiling meta-analysis of in vitro matured bovine oocyte data and in vivo matured human oocyte data to identify classifier orthologs for different melatonin environments was performed in this study. Despite the discreet procedures manifested, some limitations in result interpretation should be greatly concerned to avoid over-assumption. As described in the methodology, the in vitro and in vivo matured oocyte samples used in this study were procured from different studies with limited information about them. For revised instance, melatonin concentrations in follicles containing in vivo matured human oocytes were not indicated, while the concentrations were clearly stated in bovine oocytes' culture systems. Undetermined differences in melatonin concentrations between such environments should reduce the sensitivity of differential analysis in this study. With more available in vivo matured oocyte profiling data with well-defined melatonin environments in the future, the list of classifier orthologs could be improved for better generalization.

As expected, our previously reported procedures [START_REF] Chokeshaiusaha | Construction of pooled oocyte expression profiles of rhesus monkey and mouse for concurrent meta-analyses[END_REF] were also applicable to the selection of candidate orthologs differentially expressed between low and high melatonin environments (Fig 2 andTable 2).

While these orthologs could be used to annotate essential melatonin functions on in vivo and in vitro matured oocytes, this study aimed to identify the most reliable classifier orthologs among them. Due to such delegation, we further verified the candidate orthologs with the zero-inflated regression test to reduce the mean-variance bias (Fig 3). Encoded products of the verified candidate orthologs were associated with in various oocyte metabolisms (Table 3). In detail, ADM encoded an angiogenic factor, adrenomedullin-an important vasorelaxant required during folliculogenesis and oocyte maturation processes [START_REF] Malamitsi-Puchner | Concentrations of angiogenic factors in follicular fluid and oocyte-cumulus complex culture medium from women undergoing in vitro fertilization: Association with oocyte maturity and fertilization[END_REF][START_REF] Manau | Follicular fluid concentrations of adrenomedullin, vascular endothelial growth factor and nitric oxide in IVF cycles: Relationship to ovarian response[END_REF]. CYB561D2 encoded a member of trans-membrane (TM) protein family involving stress defense and iron metabolism [START_REF] Asard | Cytochromes b561: Ascorbate-mediated trans-membrane electron transport[END_REF]. ENDOV encoded Endoribonuclease V-an enzyme involving with the DNA repair process of oocyte [START_REF] Cao | Endonuclease V: An unusual enzyme for repair of DNA deamination[END_REF]. KIAA0586 encoded a conserved centrosomal protein required for sonic hedgehog/SHH signaling during the oocyte developmental process [START_REF] Lee | Sonic hedgehog signaling mediates resveratrol to improve maturation of pig oocytes in vitro and subsequent preimplantation embryo development[END_REF]. NR1H3-activating ligands helped with oocyte meiosis resumption [START_REF] Steffensen | Reduced fertility and inability of oocytes to resume meiosis in mice deficient of the Lxr genes[END_REF].

PLEKHM1 could regulate the autophagy process which was crucial in the regulation of oocyte metabolism [START_REF] Sagrillo-Fagundes | Melatonin: The smart molecule that differentially modulates autophagy in tumor and normal placental cells[END_REF][START_REF] Lin | Role of autophagy in modulating postmaturation aging of mouse oocytes article[END_REF][START_REF] Shen | Autophagy is required for proper meiosis of porcine oocytes maturing in vitro[END_REF]. SLC47A1 encoded an extrusion protein for drug and toxic metabolites [START_REF] Yonezawa | Importance of the multidrug and toxin extrusion MATE/SLC47A family to pharmacokinetics, pharmacodynamics/toxicodynamics and pharmacogenomics[END_REF]. SYT11 encoded synaptotagmin-11 for exocytosis and endocytosis regulation [START_REF] Wang | Synaptotagmin-11 inhibits clathrin-mediated and bulk endocytosis[END_REF]. RFT1 encoded a major folate transporter [START_REF] Umapathy | Cloning and functional characterization of the proton-coupled electrogenic folate transporter and analysis of its expression in retinal cell types[END_REF], which helped to accumulate folate for normal oogenesis and early embryogenesis [START_REF] Meredith | Growing Mouse Oocytes Transiently Activate Folate Transport via Folate Receptors As They Approach Full Size1[END_REF]. Lastly, MRPS18C encoded a mitochondrial ribosomal protein and thus was required for other mitochondrial protein syntheses [START_REF] Van Blerkom | Mitochondrial function in the human oocyte and embryo and their role in developmental competence[END_REF].

The high melatonin environment was shown to enhance expressions of all verified candidate orthologs (Fig 3 andTable 3). Supporting this notice, enhanced expression of KIAA0586 by melatonin was already notified in the other study [START_REF] Lee | Melatonin influences the sonic hedgehog signaling pathway in porcine cumulus oocyte complexes[END_REF]. Since functions of these orthologs were strongly associated with qualified oocyte development as previously mentioned, their diminished expressions in the low melatonin environment should thus adversely affect the oocyte development (Table 3). As far as we know, shreds of evidence of enhanced functions of most verified orthologs by melatonin were still not warranted in mammalian oocytes-by which future studies to confirm such effects were encouraged.

revised

In the current study, the classifying potential of the verified orthologs was determined by SVM to assess their reliability and reproducibility. According to the result, low misclassification scores and high AUCs manifested among tested learning datasets (Fig 4) hereby indicated such orthologs as potential classifiers for melatonin environments. Though utilization of cross-species data in this study partially suggested such classifier orthologs to be universally applicable among mammalian oocytes, further confirmation in other mammals with known melatonin environments were still required. Since overall expression levels of acquired classifier orthologs could imply sufficiency of melatonin, we, on a certain level proposed them as novel parameters to provide some clues about beneficial effects of melatonin supplement on target oocytes in any culture environments of interest. Such knowledge greatly relieved burdens for researchers to screen for universal melatonin effects on all ortholog expressions which were not possible by the conventional approach. It should however be noted that these classifier orthologs might not only be specific to the melatonin environment and thus required further classification tests with other possible oocyte environments.

Conclusion

Novel classifier orthologs of oocytes matured under low and high melatonin environments were identified by cross-species transcription profiling meta-analysis in this study. Such an outcome not only offered us future melatonin research topics to evaluate classifying potentials of the verified orthologs at biological levels but also the complete process to acquire them from available RNA-Seq data. We would like to note that the analytical procedures demonstrated in the current study could be applied with other oocyte environments of interest-by which identification of other novel parameters of other oocyte culture environments was possible. Predicted probability plots and ROC curves. Predicted probabilities plots were displayed by colored dots, where red dots corresponded to class 0-low melatonin environment, and green dots for class 1-high melatonin environment. The misclassification rate acquired from each learning dataset was shown under the plot (5-fold cross-validation, Bootstrap, and Monte-Carlo-cross-validation). Adjacent to each predicted probability plot, the ROC curve generated from its same learning dataset was drawn with AUC value demonstrated.

revised 1 .

 1 Introductionin vitro embryo production (IVP) is regarded as a well-established assisted reproductive technology

  . The silhouette width values acquired from PAM clustering were illustrated by the bar plot drawn by the 'ggplot2' packages. A bubble plot manifesting -log 10 FDR values and log 2 fold-change values of the verified candidate orthologs was drawn using the 'ggplot2' package. Probability plots and Receiver operating characteristic (ROC) curves were drawn to illustrate the classification accuracy of the verified candidate orthologs by the 'ggplot2' package.

  1), oocyte samples were correctly categorized into 2 clusters corresponding to their melatonin environments regardless of the species (Fig 2).

revisedFigure 1 .

 1 Figure 1. Silhouette plot of PAM clustering. The average silhouette width was 0.395 when dividing oocyte samples into 2 clusters.

Figure 2 .

 2 Figure 2. Heatmap of PAM clustering results. Rows of heatmap represented scaled expression values of the candidate orthologs clustered by k-means clustering, while the columns represented oocyte samples clustered by PAM.Figure 3. Bubble plot of the verified candidate orthologs. The blue and red colors indicated scaled expression levels of target orthologs expressed by oocytes in the low melatonin environment comparing to those obtained from the high melatonin environment. The color intensities were corresponded to the -log 10 FDR values, while the dot sizes were proportional to the absolute log 2 fold-change values of differentially expressed orthologs.Figure 4.Predicted probability plots and ROC curves. Predicted probabilities plots were displayed by colored dots, where red dots corresponded to class 0-low melatonin environment, and green dots for class 1-high melatonin environment. The misclassification rate acquired from each learning dataset was shown under the plot (5-fold cross-validation, Bootstrap, and Monte-Carlo-cross-validation). Adjacent to each predicted probability plot, the ROC curve generated from its same learning dataset was drawn with AUC value demonstrated.

Figure 3 .

 3 Figure 2. Heatmap of PAM clustering results. Rows of heatmap represented scaled expression values of the candidate orthologs clustered by k-means clustering, while the columns represented oocyte samples clustered by PAM.Figure 3. Bubble plot of the verified candidate orthologs. The blue and red colors indicated scaled expression levels of target orthologs expressed by oocytes in the low melatonin environment comparing to those obtained from the high melatonin environment. The color intensities were corresponded to the -log 10 FDR values, while the dot sizes were proportional to the absolute log 2 fold-change values of differentially expressed orthologs.Figure 4.Predicted probability plots and ROC curves. Predicted probabilities plots were displayed by colored dots, where red dots corresponded to class 0-low melatonin environment, and green dots for class 1-high melatonin environment. The misclassification rate acquired from each learning dataset was shown under the plot (5-fold cross-validation, Bootstrap, and Monte-Carlo-cross-validation). Adjacent to each predicted probability plot, the ROC curve generated from its same learning dataset was drawn with AUC value demonstrated.

Figure 4 .

 4 Figure 2. Heatmap of PAM clustering results. Rows of heatmap represented scaled expression values of the candidate orthologs clustered by k-means clustering, while the columns represented oocyte samples clustered by PAM.Figure 3. Bubble plot of the verified candidate orthologs. The blue and red colors indicated scaled expression levels of target orthologs expressed by oocytes in the low melatonin environment comparing to those obtained from the high melatonin environment. The color intensities were corresponded to the -log 10 FDR values, while the dot sizes were proportional to the absolute log 2 fold-change values of differentially expressed orthologs.Figure 4.Predicted probability plots and ROC curves. Predicted probabilities plots were displayed by colored dots, where red dots corresponded to class 0-low melatonin environment, and green dots for class 1-high melatonin environment. The misclassification rate acquired from each learning dataset was shown under the plot (5-fold cross-validation, Bootstrap, and Monte-Carlo-cross-validation). Adjacent to each predicted probability plot, the ROC curve generated from its same learning dataset was drawn with AUC value demonstrated.

Acknowledgments

The current study was fully supported by the Thailand Research Fund (TRF) through New Research Scholar Program (Grant No. TRG5880003). We would like to show our gratitude to the Office of Academic Resource and Information Technology (ARiT), Rajamangala University of Technology Tawan-OK (RMUTTO) for the contribution of the server computer for the burden analysis.

Research funding: Thailand Research Fund (TRF) through New Research Scholar Programme (Grant No. TRG5880003)

Conflict of interest

We certify that there is no conflict of interest with any financial organization regarding the material and methods discussed in the manuscript.