F Cl

Lapatinib Arpan Sasmal

Thierry Roisnel

Jitendra K Bera

Henri Doucet
email: henri.doucet@univ-rennes1.fr

Jean-François Soulé
email: jean-francois.soule@univ-rennes1.fr

Reactivity of 3-Bromofuran in Pd-Catalyzed C-H Bond Arylation Toward the Synthesis of 2,3,5-Triaryl Furans

Keywords: Palladium, Furans, C-H Bond Arylation, Aryl Bromides

The C-H bond arylation at C2 position of 3-bromofuran was performed using aryl bromides as coupling partners in the presence of phosphine free Pd(OAc) 2 /KOAc in DMA. C2,C5-di-or even C2,C4,C5-tri-arylated 3bromofuran derivatives were obtained using a larger amount of aryl bromides. Then, C2,C3,C5-triarylated furans -containing three different aryl groups-were synthesized via C2-H bond arylation, Suzuki reaction, and C5-H bond arylation sequence.

Figure 1. Relevant Drugs Containing a 2-Arylfuran Motif

Therefore, the development of efficient and regioselective synthetic pathways of multi-substituted furans has been among important research subjects in heterocyclic chemistry. Beside the traditional approaches involving cyclocondensation, cycloisomerization or multi-steps synthesis, 6 the direct functionalization of furan rings via palladium catalyzed C-H bond arylation is a modular approach to access to polyarylated furans. 7 The first example of Pd-catalyzed C2-H bond arylation of furan with activated aryl bromides was reported in 1990 by Ohta and co-workers using 5 mol% Pd(PPh3)4 associated with KOAc in DMA at 150 ºC (Figure 2a). 8 However, low to moderate yields in favor of the C2-arylated furans were often obtained due to the high reactivity of C5-H bond of C2-arylated furans leading to the formation of C2,C5-diarylated furans. 9 On the other hand, the reactivity of C3-substituted furans attracted less attention. The regioselectivity control of the C-H bond arylation of ethyl 3-furoate depends on the reaction conditions. In 2003, Sharp and co-workers reported the C2 arylation with electron-deficient aryl bromides using 5 mol% Pd(PPh3)4 associated with KOAc in non-polar solvents such as toluene at 110 °C (Figure 2b, right). 10 In contrast, they reported two examples of C5-arylation of ethyl 3-furoate using polar solvents, such as NMP, with Pd/C as catalyst (Figure 2b, left). 10 Later, Doucet and co-workers also investigated the regioselectivity of C-H bond arylation of ethyl 3-furoate and observed that using 0.2 mol% Pd(OAc)2/dppb in the presence of KOAc in DMA, C5-arylated furans were obtained as major products along with the formation of C2-arylated furans as minor products (Figure 2c). 11 Doucet and co-worker also reported one example of C2-arylation of 3-formylfuran using a similar catalytic system (Figure 2d). 12 Conversely, the chemoselective C-H bond arylations of thiophenes or pyrazoles bromo-substituted at C3 or C4 positions with aryl halides allowed the construction of polyarylated heteroaromatics. 13 To the best of our knowledge, the reactivity of 3-bromofuran in direct arylation was not reported yet. Herein, we investigate on the coupling of 3-bromofuran with electron-deficient aryl bromides and on their polyarylation to construct triarylated furans (Figure 2e). We selected 4-bromoacetophenone as the coupling partner to investigate the chemoselectivity of the coupling with 3-bromofuran. During the reaction development, we determined that the best conditions to regioselectively obtain C2-arylated 3-bromofuran 1 in 47% yield were 1 mol% of Pd(OAc)2 as catalyst, 2 equivalents of KOAc as base in DMA at 120 ºC over 8 h (

74%

OMe MeO arylation at C2 position was confirmed by X-ray-analysis of 1. 14,[START_REF]For the synthesis of 2-aryl-3-bromofurans via cyclization reactions, see: Obrecht, D[END_REF] The amount of 3-bromofuran (5 equivalents) was critical to isolate 1 in a good yield due to oligomerization as major side reaction (Table 1, entries 1 and 2). The reaction time was also a very important factor for this coupling, as shorter (4 h) or longer reaction times (16 h) resulted in lower yields (Table 1, entries 3 and 4). A higher temperature of 150 ºC led to oligomerization, while at 80 ºC a partial conversion of 4-bromoacetophenone was observed (Table 1, entries 5 and 6). Other palladium sources than Pd(OAc)2 as well as different catalyst loadings gave lower yields in 1 (Table 1, entries 7-11). Reactions performed in other solvents such as cyclopentyl methyl ether (CPME), 1,4-dioxane or pentan-1-ol failed to afford the arylated furan 1 (Table 1, entries 12-14). The concentration plays a major role to avoid the oligomerization. Indeed, 0.125 molar was determined to be the best choice, as a concentration of 0.25 molar afforded more oligomers and when the reaction was performed at 0.1 molar a partial conversion was obtained (Table 1, entries 16 and 17). Finally, our attempts to inhibit the degradation of 3-bromofuran by addition of TEMPO or catechol were unsuccessful (Table 1, entries 17 and 18). (5 equiv.) (0.5 mmol) bromides (Scheme 1). Firstly, the reactivity of a set of para-substituted aryl bromides was investigated. The reaction proceeded smoothly using aryl bromides bearing an electron-withdrawing group such as cyano, trifluoromethyl, propionyl, or chloro affording the C2-arylated 3-bromofurans 2-5 in moderate yields. Notably, when non-activated aryl bromides were employed (e.g., bromobenzene or 4bromoanisole), the corresponding C2-arylated 3-bromofurans were not obtained due to the oligomerization of 3-bromofuran. Meta-substituted aryl bromides, such as 3-bromobenzonitrile and 3bromobenzotrifluoride were chemoselectively coupled with 3-bromofuran to give the arylated compounds 6 and 7 in 50% yields, again without formation of diarylated products. The reaction was slightly sensitive to steric hindrance, as from 2-bromobenzonitrile the coupling product 8 was isolated in only 42% yield. From less bulky and activated 4-bromo-2-fluorobenzaldehyde, Pd-catalyzed C2 arylation lead to the furan derivative 9 in 52% yield. 2-Bromonaphtalene displayed a higher reactivity than 1bromonaphtalene in Pd-catalyzed C2 arylation of 3-bromofuran, owing steric factors. 3-Bromopyridine was also successfully coupled with 3-bromofuran to afford 12 in good yield. Other N-heteroaryl bromides such as 5-bromopyrimidine and 3-bromoquinoline gave the C2-heteroarylated 3-bromofuran derivatives 13 and 14 in 45% and 50% yield, respectively. Scheme 1. Scope of (Hetero)Aryl Bromides in Pd-Catalyzed Direct C2-Arylation of 3-Bromofuran C-H bond arylation of furan using palladium catalysis often occurred at C2 and C5 positions. 8 We then decided to employ an excess amount of aryl bromides to have access to C2,C5-diarylated 3-bromofurans in one step (Scheme 2). Reaction from 2 equivalents of 4-bromobenzonitrile and 3-bromofuran in the presence of 2 mol% Pd(OAc)2 associated with 3 equivalents of KOAc did not afford the desired diarylated product 15, but only mono-arylated furan 2 albeit in a poor yield. Surprisingly, in the presence of 2 equivalents of 4-bromobenzotrifluoride, 3-bromofuran was diarylated at C2 and C5 positions to give the furan 16 in 56% yield. This reactivity was also observed with 2-or 3-bromobenzotrifluorides and 2bromo-6-(trifluoromethyl)pyridine, as the diarylated furans 17-19 were selectively obtained in 47-53% yields. It is noteworthy to mention that in all cases, the formation of mono-arylated 3-bromofurans was also observed in 5-10% yields. The reason of the outperformance of aryl bromides bearing a trifluoromethyl group for the multi C-H bond arylations remains unclear. However, similar phenomena have been previously observed in the Pd-catalyzed multiple C-H bond arylations of other heterocycles (e.g., thiophenes, pyrroles, imidazoles and benzofuran). 16 Scheme 2. Scope of Hetero)Aryl Bromides in Pd-Catalyzed Direct C2,C5-Diarylation of 3-Bromofuran

We also performed the multiple C-H bond arylation of 3-bromofuran with 1-bromo-3,5bis(trifluoromethyl)benzene using 2 mol% of Pd(OAc)2 as catalyst in DMA (Scheme 3). In the presence of 2 equivalents of this aryl bromide with 3 equivalents of base, the formation of triarylated furan 20 as well as the formation of a diarylated furan was observed. Notably, we were not able to separate them by flash chromatography. Furthermore, we set-up the reaction using 4 equivalents of 1-bromo-3,5bis(trifluoromethyl)benzene to get selectively the triarylated furan 20 in good 61% yield. We also tried to prepare 2,5-diarylated 3-bromofurans with two different aryl units by iterative C-H bond arylations (Scheme 4). However, our attempts to selectively arylate the C5-H bond of 3-bromo-2-(naphthalen-1-yl)furan (10) with 4-bromobenzonitrile remained unsuccessful even in the presence of the highly active diphosphine-palladium complex [PdCl(C3H5)(dppb)]. Instead, a complex mixture of side products was observed by GC-MS analysis, which included debrominated products. Therefore, we decided to transform the C-Br bond into a less reactive C-Ph bond to later investigate the reactivity of the C5-H bond. From a mixture of 3-bromo-2-(naphthalen-1-yl)furan ( 10) and phenylboronic acid (1.5 equivalents) in the presence of 2 mol% of a diphosphine-palladium catalyst and 2 equivalents of K3PO4 in dioxane at 100 ºC over 16h, the 2,3-diarylfuran 22 was obtained in 78% yield (Scheme 5).

Scheme 5. Reactivity of C-Br Bond of 2-Naphtyl-3-Bromofuran in Pd-Catalyzed Suzuki Reaction

Only few methods are reported to prepare 2,3,5-triaryl furans containing three different aryl groups and they involve cyclization reactions. [START_REF] Dudnik | For selected examples of synthesis of 2,3,5-triaryl furans containing three different aryl groups[END_REF] To the best of our knowledge, there is no example of direct arylation of 2,3-diarylfurans. Therefore, we further demonstrated the potential of this methodology with the introduction of a third aryl group, by the sequential synthesis of 2,3,5-triarylfuran derivatives containing three different aryl units (Scheme 6). From a mixture of 2-naphtyl-3-phenylfuran (22) and 4bromobenzonitrile (1.5 equiv.) in the presence of PdCl(C3H5)(dppb) (2 mol%) associated with KOAc as base in DMA at 150 ºC, the arylated product 23 was obtained in 71% yield. The arylation took place at the most acidic C-H bond, namely C5 position. Similarly, the reaction with 3-bromoquinoline afforded the 2,3,5-triaryl furan 24 in 81% yield. Scheme 6. Reactivity of 2-Naphtyl-3-phenylfuran in Pd-Catalyzed Direct C5-Arylation

In summary, we have demonstrated that 3-bromofuran is an attractive platform for the synthesis of polyarylated furans. Firstly, conditions for the regioselective C2-H bond arylation of 3-bromofuran with activated aryl bromides have been developed using phosphine-free Pd(OAc)2 catalyst. Interestingly, the reaction was very chemoselective. There was no formation of products resulting from the activation of the C-Br bond on the furan ring and this C-Br bond remained untouched after the C2-arylation if the aryl bromide bears an electron-withdrawing group. The use of a larger amount of aryl bromide couplings partners allowed the access to C2,C5-diarylfurans and even a C2,C4,C5-triarylfuran in the case of 3,5bis(trifluoromethyl)bromobenzene. Moreover, C2,C3,C5-triarylfurans containing three different aryl units have been prepared in good yields via C2-H bond arylation, Suzuki cross-coupling reaction, and C5-H bond arylation sequence.

Experimental Section

All reactions were carried out under argon atmosphere with standard Schlenk techniques. DMA was purchased from Acros Organics and was not purified before use. 1 H and 13 C NMR spectra are recorded on Bruker AV III 400 MHz NMR spectrometer equipped with BBFO probehead. Chemical shifts (δ) are reported in parts per million relative to residual chloroform (7.28 ppm for1 H; 77.23 ppm for 13 C), constants are reported in Hertz. 1 H NMR assignment abbreviations are the following: singlet (s), doublet (d), triplet (t), quartet (q), doublet of doublets (dd), doublet of triplets (dt), and multiplet (m). All reagents were weighed and handled in air.

Procedures

General procedure for the C2-arylation reaction of 3-bromofuran: As a typical experiment, the reaction of the aryl bromide (0.5 mmol), 3-bromofuran (2.5 mmol, 5 equiv., 367 mg) and KOAc (1 mmol, 2 equiv., 98 mg) at 120 °C during 8 h in DMA (2 mL) in the presence of Pd(OAc)2 (1 mol%, 0.005 mmol, 1.12 mg,) (see tables or schemes) under argon affords the arylation product after evaporation of the solvent and purification on silica gel.

1-(4-(3-Bromofuran-2-yl)phenyl)ethan-1-one 1 :

From 4-bromoacetophenone (100 mg, 0.5 mmol) and 3-bromofuran (2.5 mmol, 5 equiv., 367 mg), 1 was obtained in 47% (62 mg) yield.

4-(3-Bromofuran-2-yl)benzonitrile 2:

From 4-bromobenzonitrile (91 mg, 0.5 mmol) and 3-bromofuran (367 mg, 2.5 mmol, 5 equiv.), 2 was obtained in 51% (63 mg) yield.

3-Bromo-2-(4-(trifluoromethyl)phenyl)furan 3:

From 4-bromobenzotrifluoride (113 mg, 0.5 mmol) and 3-bromofuran (367 mg, 2.5 mmol, 5 equiv.), 3 was obtained in 45% (65 mg) yield. Elemental analysis: calcd (%) for C11H6BrF3O (291.07): C 45.39, H 2.08; found: C 45.51, H 2.34.

1-(4-(3-Bromofuran-2-yl)phenyl)propan-1-one 4:

From 4-bromopropiophenone (107 mg, 0.5 mmol) and 3-bromofuran (367 mg, 2.5 mmol, 5 equiv.), 4 was obtained in 42% (59 mg) yield.

3-Bromo-1-(naphthalen-2-yl)furan 10: From 1-bromonaphtalene (103 mg, 0.5 mmol) and 3bromofuran (367 mg, 2.5 mmol, 5 equiv.), 10 was obtained in 38% (52 mg) yield.

5-(3-Bromofuran-2-yl)pyrimidine 13:

From 5-bromopyrimidine (79 mg, 0.5 mmol) and 3-bromofuran (367 mg, 2.5 mmol, 5 equiv.), 13 was obtained in 52% (58 mg) yield. Elemental analysis: calcd (%) for C8H5BrN2O (225.05): C 42.70, H 2.24; found: C 42.96, H 2.37.

3-(3-Bromofuran-2-yl

)quinoline 14: From 3-bromoquinoline (104 mg, 0.5 mmol) and 3-bromofuran (367 mg, 2.5 mmol, 5 equiv.), 14 was obtained in 50% (69 mg) yield. General procedure for the C2,C5-diarylation reaction of 3-bromofuran: As a typical experiment, the reaction of the aryl bromide (1 mmol), 3-bromofuran (0.5 mmol, 1 equiv., 73 mg) and KOAc (1.5 mmol, 3 equiv., 147 mg) at 120 °C during 8 h in DMA (2 mL) in the presence of Pd(OAc)2 (2 mol%, 0.01 mmol, 2.24 mg,) (see tables or schemes) under argon affords the diarylation product after evaporation of the solvent and purification on silica gel.

3-Bromo-2,5-bis(4-(trifluoromethyl)phenyl)furan 16: From 4-bromobenzotrifluoride (223 mg, 1 mmol, 2 equiv) and 3-bromofuran (73 mg, 0.5 mmol, 1 equiv.), 16 was obtained in 56% (122 mg) yield.

A.S. acknowledges the Indo-French Centre for the Promotion of Advanced Research (IFCPAR, Nº 5705-1) for a PhD grant and research support.

Figure 2 .

 2 Figure 2. Reactivity of (Substituted) Furan Derivatives in Pd-Catalyzed C-H Bond Arylation

[ a ] 1 X-ray of 1

 a11 determined by GC-MS analysis using n-dodecane as internal standard. Isolated yield is shown in parenthesesHaving determined the best conditions to selectively arylate the C-H bond at C2 position of 3-bromofuran without cleavage of the C-Br bond on the furan ring, we turned our attention to the scope of aryl

Scheme 3 .

 3 Scheme 3. Pd-catalyzed Direct C2,C4,C5-Triarylation of 3-Bromofuran with 1-Bromo-3,5bis(trifluoromethyl)benzene

Scheme 4 .

 4 Scheme 4. Reactivity of 2-Naphtyl-3-bromofuran in Pd-Catalyzed Direct C5-Arylation

1 H

 1 NMR (400 MHz, CDCl3, 25 °C): δ (ppm) = 8.09 (d, J = 8.4 Hz, 2H), 7.69 (d, J = 8.4 Hz, 2H), 7.49 (d, J = 1.7 Hz, 1H), 6.59 (d, J = 1.7 Hz, 1H).

1 H

 1 NMR (400 MHz, CDCl3, 25 °C): δ (ppm) = 8.12 (d, J = 8.4 Hz, 2H), 7.70 (d, J = 8.4 Hz, 2H), 7.50 (d, J = 1.9 Hz, 1H), 6.61 (d, J = 1.9 Hz, 1H).13 C NMR (100 MHz, CDCl3, 25 °C): δ (ppm) = 147.5, 142.6, 132.9, 128.4, 125.5 (q, J = 4.1 Hz), 125.4, 124.0 (q, J = 274.2 Hz), 116.6, 98.0.19 F NMR (376 MHz, CDCl3, 25 °C): δ (ppm) = -62.7 (s).

1 H

 1 NMR (400 MHz, CDCl3, 25 °C): δ (ppm) = 9.34 (s, 2H), 9.18 (s, 1H), 7.53 (d, J = 2.0 Hz, 1H), 6.62 (d, J = 2.0 Hz, 1H).13 C NMR (100 MHz, CDCl3, 25 °C): δ (ppm) = 157.3, 153.0, 144.0, 143.9, 124.6, 116.8, 99.5.

Table

Pd-catalyzed C2-arylation of Furan (Ohta et al.)

  1, entry 1). The regioselectivity of the

	O a. O Br Ar + Pd(PPh 3 ) 4 (5 mol%) KOAc, DMA, 150 ºC b.	O Pd(PPh 3 ) 4 (5 mol%) KOAc	Ar	Ar	O
	+	toluene, 110 ºC
	Br Ar

Pd-catalyzed C2-or C5-arylation of ethyl 3-furoate (Sharp et al.)

		EtO 2 C	EtO 2 C
		Pd(OAc) 2 /dppb
		(0.2 mol%)
	+	Br Ar
	O		O

c. Pd-catalyzed C5-arylation of ethyl 3-furoate (Doucet et al.)

					low to moderate yields	(7 examples)
	EtO 2 C			Pd/C (5 mol%)
	O	Ar		KOAc NMP, 110 ºC	(3 equiv.)
	(2 examples)		low yields		moderate yields	(5 examples)
	EtO 2 C					EtO 2 C	EtO 2 C
							+
	(2 equiv.)			KOAc, DMA, 130 ºC moderate yields	Ar	O minor	major	Ar
							(2 examples)
				1/2[PdCl(C 3 H 5 )] 2 /
					Tedicyp
		+	Br Ar	(1 mol%)
	O						O
	e. Reactivity of 3-bromofuran in Pd-catalyzed C-H bond arylation (this work)
	Br				[Pd]	Br	Br
			+	Br Ar			or
	O						Ar	O	O	Ar

d. Pd-catalyzed C2-arylation of 3-formylfuran (Doucet et al.)

	OHC	OHC	OHC
			+
	NaOAc, DMA, 150 ºC	O
	(2 equiv.)	90	10
		(1 example)

Table 1 .

 1 Reaction Development

	[a]

Entry

Variation from standard conditions Yield in 1 (%)

  H NMR (400 MHz, CDCl3, 25 °C): δ (ppm) = 7.98 -7.90 (m, 3H), 7.73 (dd, J = 1.3, 7.1 Hz, 1H), 7.61 (d, J = 2.0 Hz, 1H), 7.59 -7.53 (m, 3H), 6.68 (d, J = 2.0 Hz, 1H).

	13 C NMR (100 MHz, CDCl3, 25 °C): δ (ppm) = 150.5, 142.8, 133.9, 131.7, 130.0, 129.0, 128.5, 126.8, 126.3,
	126.1, 125.1, 115.3, 99.1.
	Elemental analysis: calcd (%) for C14H9BrO (273.13): C 61.57, H 3.32; found: C 61.46, H 3.28.
	3-Bromo-2-(naphthalen-2-yl)furan 11: From 2-bromonaphtalene (103 mg, 0.5 mmol) and 3-
	bromofuran (367 mg, 2.5 mmol, 5 equiv.), 11 was obtained in 45% (61 mg) yield.
	1 H NMR (400 MHz, CDCl3, 25 °C): δ (ppm) = 8.45 (s, 1H), 8.08 (dd, J = 1.8, 8.7 Hz, 1H), 7.94 -7.86 (m, 2H),
	7.86 -7.81 (m, 1H), 7.54 -7.45 (m, 3H), 6.59 (d, J = 2.0 Hz, 1H).
	13 C NMR (100 MHz, CDCl3, 25 °C): δ (ppm) = 149.2, 142.1, 133.3, 132.9, 128.6, 128.3, 127.8, 127.3, 126.6,
	126.6, 124.9, 123.4, 116.6, 96.6.
	Elemental analysis: calcd (%) for C14H9BrO (273.13): C 61.57, H 3.32; found: C 61.79, H 3.51.
	3-(3-Bromofuran-2-yl)pyridine 12: From 3-bromopyridine (79 mg, 0.5 mmol) and 3-bromofuran (367
	mg, 2.5 mmol, 5 equiv.), 12 was obtained in 45% (50 mg) yield.

1 1 H NMR (400 MHz, CDCl3, 25 °C): δ (ppm) = 9.23 (brs, 1H), 8.57 (brs, 1H), 8.38 -8.11 (m, 1H), 7.47 (d, J = 2.0 Hz, 1H), 7.36 (dd, J = 4.8, 8.2 Hz, 1H), 6.57 (d, J = 2.0 Hz, 1H). 13 C NMR (100 MHz, CDCl3, 25 °C): δ (ppm) = 148.8, 146.8, 146.5, 142.9, 132.5, 123.5, 116.6, 97.8. Elemental analysis: calcd (%) for C9H6BrNO (224.06): C 48.25, H 2.70; found: C 48.21, H 2.83.

  NMR (400 MHz, CDCl3, 25 °C): δ (ppm) = 9.50 (d, J = 2.3 Hz, 1H), 8.69 (d, J = 2.3 Hz, 1H), 8.11 (d, J = 8.4 Hz, 1H), 7.87 (d, J = 8.1 Hz, 1H), 7.73 (t, J = 6.8 Hz, 1H), 7.57 (t, J = 7.9 Hz, 1H), 7.52 (d, J = 2.0 Hz, 1H), 6.61 (d, J = 2.0 Hz, 1H).13 C NMR (100 MHz, CDCl3, 25 °C): δ (ppm) = 147.7, 147.2, 146.9, 143.0, 131.8, 130.0, 129.4, 128.3, 127.6, 127.4, 123.3, 116.7, 98.0. Elemental analysis: calcd (%) for C13H8BrNO (274.12): C 56.96, H 2.94; found: C 57.08, H 2.69.

1 H

H NMR (400 MHz, CDCl3,

°C): δ (ppm) = 8.08 (d, J = 8.8 Hz, 2H), 8.02 (d, J = 8.9 Hz, 2H), 7.47 (d, J = 1.9 Hz, 1H), 6.58 (d, J = 1.9 Hz, 1H), 3.02 (q, J = 7.2 Hz, 2H), 1.25 (t, J = 7.2 Hz, 3H).

Acknowledgment

We thank CNRS and ''Rennes Metropole" for providing financial support.

Supporting Information

Supporting information contains the 1 H and 13 C charts of all compounds.


Hussein Kanso

Bénédicte Quilot-Turion

Mohamed-Mahmoud Memah

Olivier Bernard

Jean-Luc Gouzé

Valentina Baldazzi

Reducing a model of sugar metabolism in peach to catch dierent patterns among genotypes

Keywords: model reduction, sensitivity analysis, structural simplication, quasi-steady-state, peach fruit, kinetic model, model calibration, gene-to-phenotype

 to simulate the accumulation of dierent sugars during peach fruit development. Two major drawbacks of this model are (a) the number of parameters to calibrate and (b) its integration time that can be long due to non-linearity and time-dependent input functions. Together, these issues hamper the use

of the model for a large panel of genotypes, for which few data are available.

In this paper, we present a model reduction scheme that explicitly addresses the specicity of genetic studies in that: i) it yields a reduced model that is adapted to the whole expected genetic diversity ii) it maintains network structure and variable identity, in order to facilitate biological interpretation. The proposed approach is based on the combination and the systematic evaluation of dierent reduction methods. Thus, we combined multivariate sensitivity analysis, structural simplication and timescale-based approaches to simplify the number and the structure of ordinary dierential equations of the model. The original and reduced models were compared based on three criteria, namely the corrected Aikake Information Criterion (AI C C ), the calibration time and the expected error of the reduced model over a progeny of virtual genotypes. The resulting reduced model not only reproduces the

Introduction

Plants are sessile organisms endowed with the capacity to alter their development, physiology, and morphology depending on the context. Plant phenotype is the result of the interaction between the environment, cultural practices and plant's genetic background (genotype). In the context of agronomy, increasing eorts have been made to select varieties that better meet consumers' expectations. Today it is clear that future breeding should account for complex plant phenotypes, responding to a large panel of criteria, including increased yield, abiotic and biotic stress tolerance, and quality of food products.

Genotype-phenotype models have been considered as the tools of the future to design new genotypes since they can help to test the performance of new genotypes (G) under dierent Environments (E) x Management (M) conditions. The challenge is to build ecophysiological models that integrate genetic information associated to specic processes (traits). In general, genotypes are dened by a set of parameters, which depends on gene expression or allelic combination, depending on the genetic complexity of the considered trait as well as the available information [START_REF] White | Gene-Based Approaches to Crop Simulation[END_REF]. Genetic-improved ecophysiological models can then be used to capture GxExM interactions. They can also be used to design ideotypes i.e. real or virtual plant cultivars expressing an ideal phenotype adapted to a particular biophysical environment, crop management, and end-use [START_REF] Letort | Parametric identication of a functionalstructural tree growth model and application to beech trees (fagus sylvatica)[END_REF][START_REF] Tardieu | Virtual plants: modelling as a tool for the genomics of tolerance to water decit[END_REF]. For this, it is necessary to combine the genetic-improved ecophysiological model with a multi-objective optimization algorithm to identify the best genotypes for specic conditions [START_REF] Quilot-Turion | Optimization of allelic combinations controlling parameters of a peach quality model[END_REF].

Construction of gene-to-phenotype models is challenging. First, the approach requires that a sole and unique model can reproduce the behavior of all genotypes, in multiple environments, the diversity observed being supported by dierent sets of parameters. Second, calibration of the models for a large number of genotypes is generally dicult, due to a large number of parameters (typically from 50 to 200 in whole-plant ecophysiological models) along with a restricted number of observations [START_REF] Martre | Chapter 14 -model-assisted phenotyping and ideotype design[END_REF][START_REF] Bertin | Under what circumstances can process-based simulation models link genotype to phenotype for complex traits? case-study of fruit and grain quality traits[END_REF]. Due to the model complexity and non-lineairities, evolutionary and bio-inspired algorithms are increasingly used both for parameter estimation and ideotype design. These methods can explore high-dimensional parameter space eciently but they rely on a large number of model evaluations, that can rapidly increase the computational time required to nd a solution. Third, the genetic architecture of complex traits can be very complex, due to epistatic and pleiotropic eects. In this sense, the presence of biologically-meaningful parameters can considerably help the interpretation of the resulting genetic architecture, facilitating the breeding process. Ideally, most the model is close to omics data, the easier the linkage between the parameters and the underlying physiological processes.

Kinetic modeling has been successfully applied to several metabolic pathways in plants [START_REF] Curien | The music industry in the digital era: Toward new contracts[END_REF][START_REF] Nägele | Mathematical modeling reveals that metabolic feedback regulation of snrk1 and hexokinase is sucient to control sugar homeostasis from energy depletion to full recovery[END_REF][START_REF] Beauvoit | Model-assisted analysis of sugar metabolism throughout tomato fruit development reveals enzyme and carrier properties in relation to vacuole expansion[END_REF]. In this spirit, a kinetic model of sugar metabolism has been developed in [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] to simulate the accumulation of dierent sugars during peach fruit development. The model correctly accounts for annual variability and the genotypic variations observed in ten genotypes derived from a larger progeny of inter-specic peach cross. At term, the objective of the research is to integrate the genetic control of sugar metabolism in this kinetic model and develop a methodology to design ideotypes by virtual breeding. To achieve this, it is necessary to estimate accurately the values of the inuential parameters of the model for the whole progeny of 106 genotypes for which few data are available. Unfortunately, the size of the parameter space and the non-linearity of the reaction rates make the calibration of the model unreliable and time-consuming.

One way to face these weaknesses is to reduce the complexity of the model [START_REF] Okino | Simplication of mathematical models of chemical reaction systems[END_REF]. Several reductions and approximation approaches exist in the literature, each one addressing a specic aspect of model complexity [START_REF] Gorban | Model reduction and coarse-graining approaches for multiscale phenomena[END_REF][START_REF] Snowden | Methods of model reduction for large-scale biological systems: a survey of current methods and trends[END_REF]. A number of methods, such the lumping method [START_REF] Wei | Lumping analysis in monomolecular reaction systems. analysis of the exactly lumpable system[END_REF][START_REF] Sunnåker | A method for zooming of nonlinear models of biochemical systems[END_REF] or the classical quasi-steadystate (QSS) approaches, aim at reducing the number of variables based on chemical or time-scale considerations [START_REF] Schauer | Quasi-steady-state approximation in the mathematical modeling of biochemical reaction networks[END_REF][START_REF] Heinrich | The regulation of cellular systems[END_REF]. Methods from sensitivity analysis may help to reduce the parameter space by identifying non-inuential parameters, whose values can be xed by broad literature data [START_REF] Turányi | Sensitivity analysis of complex kinetic systems. tools and applications[END_REF][START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF][START_REF] Vanuytrecht | Global sensitivity analysis of yield output from the water productivity model[END_REF][START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF].

Last but not least, the structure of the model itself can be simplied. Meth-ods for model decomposition [START_REF] Holme | Subnetwork hierarchies of biochemical pathways[END_REF][START_REF] Anderson | Model decomposition and reduction tools for large-scale networks in systems biology[END_REF][START_REF] Sun | Model reduction and parameter estimation of non-linear dynamical biochemical reaction networks[END_REF] aim to separate the system into sub-networks or sub-models, that are easier to analyze and parameterize.

The choice of reaction kinetics is also very important for model complexity.

In this perspective, the use of simplied enzyme kinetics [START_REF] Wang | Kinetic modeling using S-systems and lin-log approaches[END_REF][START_REF] Nikerel | Model reduction and a priori kinetic parameter identiability analysis using metabolome time series for metabolic reaction networks with linlog kinetics[END_REF][START_REF] Schmidt | Complexity reduction of biochemical rate expressions[END_REF] may be useful to avoid the emergence of numerical and identiability issues.

Dierent reduction methods can be combined together. In [START_REF] Liebermeister | Biochemical network models simplied by balanced truncation[END_REF] for instance, model decomposition is associated to variable transformation, resulting in a low-dimensional description of the exterior part of the system, whereas in [START_REF] Sunnåker | A method for zooming of nonlinear models of biochemical systems[END_REF] time scale analysis is used to identify a cluster of fast variables to be lumped together.

In the work of Apri et al. [START_REF] Apri | Complexity reduction preserving dynamical behavior of biochemical networks[END_REF] dierent reduction steps (parameter removal, node removal, variable lumping) are sequentially tested following a practical scheme: at each step, if the reduced model, after parameter reestimation, can reproduce some target outputs, the modication is selected, and rejected otherwise. From the point of view of genetic applications, a major drawback of the approach of Apri et al. [START_REF] Apri | Complexity reduction preserving dynamical behavior of biochemical networks[END_REF] is that the selection of acceptable reduction results depends on the specic target dynamics.

As a consequence, dierent target outputs (i.e. genotypes) can give rise to reduced models with dierent structures or parameters number, making their comparison dicult in the perspective of genetic studies.

The objective of this work was to provide a method to build a reduced model that is adapted to the specicity of genetic studies in that: i) it yields a reduced model that is adapted to the whole expected genetic diversity ii) it maintains network structure and variable identity, in order to facilitate the biological interpretation of the reduced model.

Similarly to the approach of Apri et al. [START_REF] Apri | Complexity reduction preserving dynamical behavior of biochemical networks[END_REF], our reduction strategy tests dierent methods in several parallel steps that, if retained, are combined together into a nal reduced model (Fig. 1).

First, multivariate sensitivity analysis was attempted to reduce the parameter space [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF]. Second, we tried to simplify the structure of the model by reducing non-linearity and time-dependent forcing, and nally, a quasisteady-state approximation based on time-scale separation was tested to reduce the size of the system. Particular attention was devoted to the systematic evaluation of the dierent reduction methods. Three main criteria were used to assess the interest of the reduction: i) the corrected AIC value, evaluating the relative gain between model simplication and loss of accuracy over an experimental dataset, ii) the calibration time, as a measure of model eciency, iii) the expected error between the original and the reduced model over a population of virtual genotypes, as a measure of the reliability of the simplication scheme.

As a case study, the proposed reduction scheme was applied to the model of sugar metabolism proposed by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF]. The resulting reduced model correctly reproduces data on the original 10 genotypes with only 9 estimated parameters (out of 14 in the original model) and a gain in calibration time over 40%. In addition, the reduced model was successfully calibrated on 30 new genotypes of the same inter-specic peach progeny, for which fewer data points were available.

The paper is organized as follows. In the next section, we briey present the original model of sugar metabolism developed by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF]. Section 3 is devoted to the description of the individual reduction methods, whereas Sections 4 and 5 present, respectively, the datasets and the numerical methods used for the assessment of the proposed model reduction. The results of the application of our reduction scheme to the model of sugar metabolism are reported in section 6. A general discussion on the advantages and limitations of our approach closes the paper.

Description of the peach sugar model

The model developed by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] describes the accumulation of four dierent sugars (sucrose, glucose, fructose, and sorbitol) in peach fruit during its development over a progeny of ten peach genotypes with contrasting sugar composition. The fruit was assumed to behave as a single big cell with two intra-cellular compartments, namely the cytosol and the vacuole. Carbon enters the fruit from the plant sap which is transformed by a metabolic network, including enzymatic reactions and transport mechanisms between the cytosol and the vacuole.

The developed dynamical model made explicit use of experimental data to describe the evolution of the sub-cellular compartment (due to fruit growth) and enzyme activities (due to fruit developmental program) over time. To this aim, measured fruit dry and fresh masses and enzyme activities were represented by genotype-specic temporal functions and provided as input to the model.

From a mathematical point of view, the model can be described as a set of parametric ordinary dierential equations:

dx dt = f (x(t), I(t), v(t), p), (1)
x(t 0 ) = x 0 ,

where t is the independent time variable in days after bloom (DAB); x ∈ R 10

is the concentration vector of metabolites in the corresponding intra-cellular compartment and x 0 ∈ R 10 in Eq.( 2) is the vector of the corresponding initial values. I ∈ R is the time-dependent input of carbon from the plant and v ∈ R

Model reduction methods

In this section, we present a reduction scheme explicitly dedicated to genetic studies that combines dierent methods in several parallel steps as shown in (Fig. 1) and explained in the next subsections.

Multivariate sensitivity analysis

Generally, in the case of complex models, estimating parameters requires a lot of eort and is known to be a dicult and challenging task. In particular, it is tricky to determine which parameters can be xed. The global sensitivity analysis methods allow to explore the inuence of each parameter on model outputs and thus to identify the key parameters that aect model performance and play important roles in model parameterization, calibration and optimization [START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF]. Multivariate sensitivity is a method developed by Lamboni et al. [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] that allows the application of global sensitivity analysis to models having a multivariate (eg. dynamic) output. The idea is to perform a principal components analysis on the outputs, and then compute the sensitivity indexes for each principal component. The results are summarized by the generalized sensitivity indices (GSI) that provide a unique ranking of the parameters over the whole output.

This method was applied to the 23 parameters of the original model and to the measured enzymatic activities v. Each parameter was studied at three levels, corresponding to 0.05, 0.5 and 0.95 quantiles of the previously estimated 14 parameters values [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] and to a variation of -20% and +20%

of the xed values for the other parameters. For time-dependent enzyme activities, the same -20% and +20% variation was applied on their average values over the whole dynamics.

In order to evaluate the impact of the genotype choice on the results of the sensitivity analysis, simulations were performed according to a factorial design, following the ANOVA model genotypes × (p 1 + . . .

p 23 + v 1 + . . . + v 7 ) 2 .
The package "Planor" in R (R Development Core Team 2015) was used.

The minimum resolution of the plan was xed by using the tool MinT [START_REF] Schürer | Mint: A database for optimal net parameters[END_REF] to test all main eects and interactions. The factorial design resulted in 10 × 3 9 = 196 830 simulations.

Multivariate sensitivity analysis was performed independently on the dynamics of the four output sugars (i.e sucrose, glucose, fructose, and sorbitol) that compose peach fruit. In order to determine the least sensitive parameters, the whole sugar phenotype has to be taken into account, with respect to the relative proportions of each sugar. For this aim, an aggregate generalized sensitivity index (aGSI ) was constructed for each parameter as

aGSI = 4 i=1 GSI i β i ( 3
)
where GSI is the generalized sensitivity indice computed for the sugar i and β i the relative proportion of sugar i in the fruit. β = (0.72, 0.13, 0.09, 0.05)

for sucrose, glucose, fructose, and sorbitol, respectively.

Structural simplication methods

This section aims to simplify the structure of the model in terms of network and reaction rates while preserving its predictive ability. The structural simplication includes the three following strategies:

Simplifying the description of enzymatic capacities

Seven enzymatic capacities V max are represented in the original model.

Some of these capacities were assumed to vary over time (temporal eect)

and/or to depend on the phenotypic group (phenotype eect), according to experimental evidences [START_REF] Desnoues | Proling sugar metabolism during fruit development in a peach progeny with dierent fructose-to-glucose ratios[END_REF]. The characteristics of enzyme capacities are summarized in Table 1. In order to simplify the model, we systematically tested the impact of the suppression of the phenotype and/or the temporal eect on each single capacity. Depending on the characteristics of the considered enzyme (Table 1), the procedure is slightly dierent:

Phenotype effect : V 1 max V 2 max → V 1 max + V 2 max 2 (4)
Temporal effect :

V max (t) → < V max (t) > t (5)
Double effect : (4) then ( 5) applied

where < . > t stands for temporal average over the whole dynamics.

u(x, t) = V max x(t) K m + x(t) (7)
where V max is the enzymatic capacity. K m is the anity of the enzyme for the substrate, x(t) is the concentration of the substrate at time t.

The objective here is to simplify Eq.( 7) in order to improve the eciency of the numerical simulation. Depending on the relative levels of the substrate concentration and the MM equation anity, two simplications of the ows' equations can be made:

Case 1: if x(t) << K m
Substrate concentration is small compared to the anity of the enzyme for the substrate then we can write: u(x, t) = Vmax Km x(t).

Case 2: if x(t) >> K m Substrate concentration exceeds the anity of the enzyme for the substrate, so that the enzyme can be supposed close to saturation: u(x, t) = V max .

Futile cycle removal

The presence of internal cycles within a metabolic network can lead to the appearance of thermodynamically unfeasible loops i.e. reactions that run simultaneously in opposite directions (for example Fig. 2) and have no overall eect on the exchange uxes of the system. This is an undesirable situation that causes numerical issues and makes the estimation of the corresponding parameter values an ill-posed problem. In this context, our strategy was to remove each futile cycle by replacing the antagonist reactions by a single eective reaction preserving the net exchange ux of the system. Dierent kinetics can be tested for the eective reaction, as alternative reduction approaches. Consistently with the previous reduction method, we decided to test two linear reaction forms, namely

u(x, t) = k i x i -k j x j (8)
and

u(x, t) = k i (x i -x j ) (9)
where x i , x j are the variables involved in the futile cycle and k i , k j are the coecients to be estimated.

Time-scale analysis and QSS approximation

Biological systems are often characterized by the presence of dierent time scales (seconds, hours, days). Following Heinrich and Schuster [START_REF] Heinrich | The regulation of cellular systems[END_REF], an appropriate measure of the time scales involved is given by

τ i (t) = - 1 Re(λ i (t)) (10)
where Re(λ i ) are real parts of the eigenvalues λ i of the Jacobian matrix of the system, along a given trajectory. The presence of fast modes in the system allows the reduction of the number of variables based on a quasi-steady-state assumption.

Based on the above information and on the analysis of time-series of the full model, variables can be divided into two groups x = (x (1) , x (2) ), where

x (1) and x (2) correspond respectively to the slow and fast variables of the system [START_REF] Heinrich | The regulation of cellular systems[END_REF][START_REF] López Zazueta | Analytical reduction of nonlinear metabolic networks accounting for dynamics in enzymatic reactions[END_REF].

Application of the QSS approximation states that dx (2) dt = f 2 (x (1) , x (2) ,

I(t), v(t), p) = 0 → x (2) ss = g(x (1) ) (11)
It follows that, after a relaxation period, the system can be approximated by the reduced model: 1) , g(x (1) ),

dx (1) dt = f 1 ((x (
I(t), v(t), p) (12)
of lower dimension.

Experimental and articial data 4.1. Experimental data

The 106 peach genotypes used in this study come from an inter-specic progeny obtained by two subsequent back-crosses between Prunus davidiana (Carr.) P1908 and Prunus persica (L.) Batsch `Summergrand' and then `Zephyr' [START_REF] Quilot | Analysis of genotypic variation in fruit esh total sugar content via an ecophysiological model applied to peach[END_REF]. They were planted in 2001 in a completely randomized design in the orchard of the INRAE Research Centre of Avignon (southern France).

Experimental monitoring of peach fruit growth and quality has been conducted in 2012, as described in [START_REF] Desnoues | Proling sugar metabolism during fruit development in a peach progeny with dierent fructose-to-glucose ratios[END_REF]. The concentration of dierent metabolites, namely sucrose, glucose, fructose, sorbitol, and hexoses phosphates, the fruit esh fresh weight and dry matter content were measured at dierent time points during fruit development, for all genotypes. In addition, the temporal evolution of enzymatic capacities (maximal activity) of the twelve enzymes involved in sugar metabolism was measured over the whole population [START_REF] Desnoues | Proling sugar metabolism during fruit development in a peach progeny with dierent fructose-to-glucose ratios[END_REF]. The resulting dynamic patterns were analyzed and compared by means of a generalized mixed linear-eect model (GLMM). Accordingly, some enzyme activities were shown to vary over time and/or depend on the phenotypic group [START_REF] Desnoues | Proling sugar metabolism during fruit development in a peach progeny with dierent fructose-to-glucose ratios[END_REF].

Training set

The 10 genotypes already used by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] were selected as the training set for our reduction strategies. They include ve genotypes having a `standard phenotype', namely a balanced fructose-to-glucose ratio at maturity between 0.6 and 0.9, and ve considered to have a `low fructose phenotype' due to the lower proportion of fructose compared with glucose based on their sugar composition at maturity [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF]. For these 10 genotypes, 3 biological measurements are available at 6 dates after bloom.

The training set was used to test each reduction method individually as well as their combination, based on the AIC C value and the calibration time (see section 5.3).

Validation set

The quality of the nal reduced model was evaluated by calibration on a validation set for which fewer data points were available (one single biological measurement at 6 dates). The idea was to select 30 additional genotypes of the inter-specic peach progeny, which in complement to the training set, represented the greatest diversity in terms of growth rate and duration. For this aim, experimentally measured growth curves were interpolated with a smoothing spline algorithm [START_REF] Chambers | Statistical models in S[END_REF] with 16.4 degrees of freedom in R (R Development Core Team 2015) and the maximum and average growth rate quantied as the maximum and the average of the growth curve's derivative over fruit development. A principal component analysis (PCA) was performed on growth rate and growth duration for the whole progeny of 106 genotypes using the R ADE4 library. The rst two principal components accounted for more than 90% of the genetic diversity. The rst axis was mainly related to the growth rate whereas the second one reected the duration of growth. As shown in Fig. 3, the ten genotypes of the original study provided a good representation of the observed diversity in growth rate. However, their growth duration was relatively short, compared to the existing variability. As a consequence, most of the new genotypes have been selected in the upper-left panel of the plan, in order to capture the greatest genetic diversity in terms of fruit development. An equal proportion of the two phenotypic groups was maintained.

Virtual genotypes

In addition to the training set, a virtual experiment was performed to evaluate the reliability of the reduction methods to variations in parameter values, initial conditions, and input functions, expected in large genetic populations. For this aim, 20 000 virtual genotypes were generated by randomly assigning model parameters and inputs, based on data from the 10 proles used in [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF].

The values of the parameters p were taken randomly using a uniform distribution between the minimum and the maximum of the previously estimated values over the set of 10 genotypes [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF]. Initial conditions, such as initial f ruit weight, and initial sugar concentration were assigned randomly using a uniform distribution within the range of observed values plus a variation of 40%.

Given the high correlation among parameters describing fruit growth curves [START_REF] Barrasso | Model-based qtl detection is sensitive to slight modications in model formulation[END_REF], model inputs, such as f ruit weight, were randomly assigned using a uniform distribution picking one of the observed growth dynamics and adding an overall random variation between zero and 10% on fruit weight. Finally, shifts in the duration of fruit development among genotypes were also considered. The maturity date was chosen randomly using a uniform distribution within the range of observed dates broaden of 40%.

Numerical methods

Parameter estimation

In this section, we aim to estimate the parameters of the models to t our observations i.e. our measured sugars concentrations. For this purpose, we note

X (k) = (X (k) 1 , . . . , X (k)
N ) the vector of the experimental observations at several times for the genotype k and suppose that:

E(X (k) i ) = M p (k) (x (k) i ) where x (k) i = (x (k) (t i )) is the set of system variables at (t i ) i∈[1,N ] , p (k)
is the vector of parameters to be estimated and M p (k) is the mathematical function relying the considered model to the data (see Appendix A for more information). Here, the observations X (k) are assumed to follow a Gaussian

law N (M p (k) (x (k) ), σ 2 k ) with constant variance σ 2 k .
The estimation of our parameters can be performed through the maximization of the likelihood. We note (p (k) , σ 2 k ) the log-likelihood function for the genotype k.

Under the assumption of observation independence, the log-likelihood can be dened as follows:

(p (k) , σ 2 k ) = - N 2 log(2π) - N 2 log(σ 2 k ) - 1 2σ 2 k N i=1 (X (k) i -M p (k) (x (k) i )) 2 (13)
A maximum log-likelihood estimator (p (k) , σk 2 ) of (p (k) , σ 2 k ) is a solution to the maximization problem:

(p (k) , σk 2 ) = arg max p (k) ,σ 2 k (p (k) , σ 2 k ) (14)
In this Gaussian case, the maximum log-likelihood estimator is thus equivalent to the ordinary least-square estimator:

p(k) = arg max p (k) N i=1 (X (k) i -M p (k) (x (k) i )) 2 (15)
σ2 k = 1 N N i=1 (X (k) i -M p(k) (x (k) i )) 2 (16)
Matlab software (MATLAB R2018a, The MathWorks Inc., Natick, MA) was used for model integration (solver ode23tb [START_REF] Hosea | Analysis and implementation of tr-bdf2[END_REF]) and calibration. A genetic algorithm (function ga [START_REF] Goldberg | Genetic algorithms in search, Optimization, and Ma-chineLearning[END_REF] of Global Optimisation Toolbox) was used for maximization of Eq. ( 15). The population size, the maximum number of generations, and the crossover probability have been respectively set at 200, 300, and 0.7. For each reduced version of the model (individual or combined reduction methods), free parameters were numerically re-estimated.

The tting process was considered at convergence when the average relative change in the best-cost function, i.e. the sum of squared errors, value over generations was less than 10 -6 . For each genotype k and reduced model, estimations procedure has been repeated ten times to take into account the stochastic nature of the genetic algorithm and to ensure the good exploration of the parameters' space. The solution having the best score was kept for subsequent analyses.

Model selection

Individual and combined reduction methods were evaluated according to three criteria of major importance for our application: the corrected Akaike Information Criterion (AI C C ), the gain in calibration time (%) and the expected error (%) between the original and reduced models.

Akaike Information Criterion

The AIC gives information on the likelihood of the proposed model based on available experimental data and weighted by the number of free parameters: [START_REF] Burnham | Model selection and multimodel inference : a practical information-theoretic approach[END_REF]:

AIC(p) = -2 (p, σ 2 ) + 2n p (17)
where n p is the number of estimated parameters p and (p, σ 2 ) is the maximum log-likelihood. In this paper, we used the corrected AIC as we deal with a small set of observations and a considerable number of parameters.

AIC C (p) = AIC(p) + 2n p (n p + 1) N -n p -1 (18)
where N is the number of observations. For genotype k and for each reduction method, we dened

∆ (k) AIC C (p (k) , p (k) ) = AIC C reduced (p (k) ) -AIC C original (p (k) ) (19)
as the AIC C dierence between the reduced and the original model. Note that ∆ AIC C is always computed using the best estimated parameter solution for the considered model. Whenever the average over the 10 genotypes (< ∆ AIC C > G ) was negative, the reduction method was validated.

Gain in calibration time

We used the calibration of a specic genotype (E43) as a proxy of the maximum expected calibration time on the population. Genotype E43 was selected because it required a long calibration time on the original model proposed by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] (approximately 11 hours on average on a 3.1GHz Intel(R) Xeon(R) processor) but it did not suer from numerical instabilities, that could complicate the calibration process. Note that the overall calibration time of a model depends both on the integration time of each evaluation step and on the convergence of the cost function that sets the actual number of generations performed by the algorithm. Both aspects may be aected by the model reduction.

To evaluate the gain in calibration time due to model reduction, parameter estimation was performed for each reduction method, following the general procedure (see section 5.2), and compared to the calibration time obtained for the original model. An initial population P 0 was randomly selected assuming a uniform distribution in the parameter range and then kept xed for all calibration processes (both original and reduced models). For models having a reduced number of parameters, the initial population was directly derived from P 0 .

The gain (G t ) was dened as the gain (in %) in calibration time T between the original and the reduced model:

G T = T original -T reduced T original × 100

Expected error

Simulations of the original and reduced models were compared by the Normalized Root Mean Square Error over the 10 model variables :

J i (p (k) , p(k) ) = 1 N S N S j=1 (x i (t j , p (k) ) -xi (t j , p(k) )) 2 max j (x i (t j , p (k) )) -min j (x i (t j , p (k) )) ∀i ∈ {1, . . . , 10} (20)
where x(t, p k ) and x(t, pk ) are the concentration predicted by the original and reduced model, respectively. Parameters for the reduced model were derived from the values of the corresponding parameters in the original model.

The quality of the QSS approximation was assessed by computing J i for each variable in the model, over the whole dynamics.

In the context of the virtual experiment, the Expected Error (%) of the reduced model was dened as the average distance J over the virtual population:

Expected Error = 1 N V G N V G k=1 < J i (p (k) , p(k) ) > ×100 (21)
with

< J i (p (k) , p(k) ) >= 1 10 10 i=1 J i (p (k) , p(k) )
where N V G is the number of virtual genotypes and 10 is the number of variables. In our case, N V G = 20 000. The Expected Error was used to quantify the reliability of the reduction.

6. Results

Strategy 1: Identication of low sensitive parameters

The objective of the sensitivity analysis was to identify parameters having a signicant inuence on the outputs of the model, over the whole dynamics and for all tested genotypes. A multivariate sensitivity analysis [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] was used for this purpose. The aggregate generalized sensitivity indices (aGSI) (see section 3.1) shown in Fig. 4 give a common ranking of model parameters according to their inuence on the whole sugar phenotype, as it is made up by the four output sugars (sucrose, sorbitol, glucose, and fructose). Parameter (p 1 ) related to the action of cell-wall invertase in fruit apoplasm and the coecient of sucrose import (p 8 ) are the most important parameters, followed by the activities of acid invertase (p 2 ), the activities of Fructokinase (p 3 ), Hexokinase (p 4 ) and the resynthesis rate of sucrose from hexose phos- phate (p 7 ). Indeed, p 1 , p 3 , and p 4 parameters are the most sensitive parameters for sucrose, fructose and glucose concentrations respectively (see Fig.

B.2 ).

Interestingly, the genotype factor is ranked third, meaning that it does not aect parameters' sensitivity as much as expected. A closer look at the Firstly, we tried to remove the temporal and the phenotype eects in the enzyme activities, v 2 , . . . , v 7 (v 1 has neither phenotype nor temporal eects).

The results of this simplication are shown in Table 2. The elimination of the phenotype eect for v 3 , v 6 and v 7 resulted in a decrease of the AIC C value for nine genotypes, neutral for one genotype, and was thus selected for the nal reduction. The elimination of the temporal eect for v 2 , v 4 , v 5 , v 7

was also advantageous on the corrected AIC results for all ten genotypes.

Nevertheless, when we tried to eliminate the temporal eect of v 7 , the resulting ∆ AIC C was positive for most genotypes. This is in line with the results of multi-variate sensitivity analysis according to which v 2 , . . . , v 6 have a low sensitivity on the four outputs of the model, whereas v 7 has a non-negligible eect on the dynamics of glucose concentration. According to these results, the elimination of the temporal eect was validated only for v 2 , v 4 , v 5 , v 6 .

In support of this choice, the test with the virtual genotypes shows that the expected error between the reduced and the original model is small (Table 2).

In the second phase, we tested the possibility of simplifying the enzymatic reaction rates (Eq.( 7)). For each reaction in the model, Fig. 5 compares the order of magnitude of the substrate x(t) to the corresponding anity K m .

The boxplots show that (Case 2, see section 3.2.2) simplication strategy can be applied only for the reaction rates u 5 and u 7 . Therefore, their reaction rates can be written as u = V max . All other ows verify the (Case 1, see section 3.2.2) and can therefore be expressed as u = Vmax Km x(t). The rates simplication improves the corrected AIC for eight genotypes and yields a substantial gain in the calibration time. The expected error over the virtual progeny is higher than in the previous reduction steps, but still in the range of accuracy of the original model [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF]. According to these observations, the enzymatic reaction rates simplication strategy was validated. and passive reversible transport (u 6 , u 8 ). Simulations showed that, whenever the genotype, the net ux mostly pointed in the direction of an export for both fructose and glucose from the vacuole to the cytosol [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF]. However, futile cycles occurred due to the presence of the active transport mechanism, that continually brings glucose and fructose back into the vacuole. Indeed, u 5 and u 6 (respectively u 7 and u 8 ) had the same evolution over the whole dynamics for all ten genotypes (Fig. 6): the active and passive transport ran simultaneously in two opposite directions.

According to our strategy (section 3.2.3), we tried to remove futile cycles by replacing reactions (u 5 , u 6 ) (respectively (u 7 , u 8 )) with an eective reaction rate of the form p 10 x 9 -p 11 x 4 (respectively p 9 x 8 -p 12 x 3 ) preserving the net export ux from vacuole to the cytosol. We compared the performance of the reduced model with respect to the original one (Table 2). The corrected AIC values were generally slightly negative, with the exception of genotypes E 1 and E 43 , suggesting an overall improvement of the model structure. Notice that the present strategy did not reduce the total parameters number but decreased model complexity and improved the calibration time.

As a further simplication, we then tried to use a special case of the above mentioned reaction rate with p 10 = p 11 (respectively p 9 = p 12 ). This time, the simplication was fully validated by the corrected AIC on all genotypes (Table 2, Eq.( 9)). The expected error over the virtual genotypes was estimated to 13% and the calibration time was lowered by 24% with respect to the original model, thanks to structural simplication and the reduction of the number of parameters to be estimated. Accordingly to these results, the simplication by Eq.( 9) was validated. 2). The expected error over a large progeny was estimated around 20%, close to the performance of the original model. On the basis of this intermediate reduced model, time scale analysis was performed to detect the possible presence of fast modes in the system. The analysis of the Jacobian matrix, indeed, conrmed the presence of dierent modes, with typical time scales spanning a few seconds up to days, for all tested genotypes (Fig. 7).

A fast transient dynamics, followed by a slow one, was observable in the numerical simulations of the original and intermediate reduced models for the hexose phosphates concentration (variable x 5 , see supplemental information, Fig. D.4). In addition, following the method proposed in [START_REF] López Zazueta | Analytical reduction of nonlinear metabolic networks accounting for dynamics in enzymatic reactions[END_REF][START_REF] Heinrich | The regulation of cellular systems[END_REF],

we analyzed the predicted concentration of sugars in both intracellular compartments, for all genotypes. The concentration of the hexose phosphate (x 5 ) was systematically lower than the concentrations of the other variables in the system, as expected for the fast components of the system (Fig. 8).

Figure 9: Normalized Root Mean Square Errors J i , i ∈ {1, . . . , 10} between the intermediate and reduced models after application of the QSSA to x 5 . The boxplot shows the variability of J i over the training set

Evaluation of the reduced model

The validity of the reduced model was veried on some new genotypes of the inter-specic peach progeny, for which few data were available. The reduced model was then calibrated on the dynamics of sugar concentration of these selected genotypes, as described in section 5.2. The results presented in (Fig. 10) showed a satisfactory agreement between model and data, all over fruit development, for most genotypes. The average N RM SE (Table D.6) ranged from 10% to 30% for the main sugars, in good agreement with estimations over the virtual progeny. These results conrmed that the reduced model oered a quality of prediction close to the original one with fewer parameters to be estimated and shorter integration time.

From a biological perspective, an important prediction of the model developed by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] was that a dierence in fructokinase anity could be at the origin of the phenotypic dierence observed between standard and low fructose genotypes.

We checked if the estimations obtained with the reduced model still supported this hypothesis. Fig. 11 shows a signicant dierence of estimated fructokinase anity between the two phenotypic groups, in agreement with the original model based on the Student t-test (p-value <2.0187e -9 )

Figure 11: Dierence in the estimated fructokinase anity between standard and low fructose phenotypes, for forty genotypes (training and validation sets). The dierence is signicant with a p-value < 2.0187e -9 .

Discussion

Models of metabolic systems are usually very complex. Complexity stems from the number of components and the high degree of non-linearity included in both the network structure and the individual reaction rates. As a consequence, metabolic models usually suer from numerical and identiability issues that seriously hamper their application in the context of genetic studies, especially when they have to be calibrated for hundreds of genotypes. In this paper, we present a reduction scheme that explicitly accounts for genomic diversity. Our approach is based on the systematic evaluation of dierent reduction methods, that, if successful, are then combined together to yield the nal reduced model. When applied to the model of sugar metabolism developed by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] our approach led to a reduced model that could be eciently calibrated on a large diversity of genotypes, for which few data are available. The reduced model showed comparable predictions and biological interpretation as the original model, with only a limited number of estimated parameters. Indeed, calibration time was reduced by 40%, a considerable improvement when considering that the calibration of the original model could span up to 30 hours for a single genotype. Moreover, mitigation of model non-linearities can help limiting numerical issues and increase the reliability of estimated parameters, an important aspect in the context of genetic studies, where large genetic populations have to be calibrated.

The proposed reduction scheme is especially suitable for dynamical models of metabolic and biochemical networks, in which a large number of chemical reactions interact with similar non-linear kinetics. In these systems, indeed, the connectivity properties of the network usually prime over the precise description of the individual rate laws [START_REF] Barabási | Network biology: understanding the cell's functional organization[END_REF]. The presence of saturating kinetic functions (like the classical Michaelis-Menten), in particular, allows the simplication of the rate function depending on the substrate range whereas the presence of redundant or opposite reactions opens the way to structural simplication of the system. The extension of these reduction steps to another kind of models is less straightforward. Crop models for instance can involve a large variety of process kinetics, one for each described physiological process. The complexity of the cellular network is replaced by the interaction of a comparatively small number of processes but described by complicated, ad-hoc kinetic functions that can involve several model components as well as external environmental variables (temperature, humidity, light). The simplication of individual rate laws is still possible but it involves case-by-case study.

Although the application of specic reduction methods is tailored to model structure, the proposed evaluation strategy is pretty generic and easily adaptable to a large range of biological models. The main objective of this work was to provide a method to build a reduced model that is adapted to the application to a large panel of genotypes. In this sense, we do not look for the best model for a given genotype but rather for the best compromise in terms of accuracy and eciency over a large genetic diversity. The question recalls the one of "model validation domain" i.e. the ability for a given model to describe data obtained in conditions dierent from those in which the model itself was calibrated [START_REF] Mairet | Twelve quick tips for designing sound dynamical models for bioprocesses[END_REF]. Here it is about selecting for a reduced model having a large validation domain and able to cope with changes in model's inputs, parameter values, and initial conditions.

For this aim, we proposed a criterium based on the simulation of a large number of virtual genotypes and the systematic comparison of the expected distance between the original and the reduced models. Virtual genotypes are built based on the variability observed in a sub-sample of the population, plus a basal variability, expressed as a random eect, to limit the bias due to the choice of the initial sample and to assure a minimal diversity across the virtual population. A few remarks are needed. First, the above method tests the reliability of the reduction, assuming that the original model is valid. In this sense, the amplitude of the basal random eect should be subject to an expert knowledge so to avoid biologically unreasonable situations, that fall outside the conditions of applicability of the model. Second, it is worth to notice that, given the virtual nature of our comparison, the reduced model is parameterized using parameter values that are directly derived from the parameters of the original model, to which it is compared. In this sense, the 'expected NRMSE error' of the reduced model represents an upper bound of its actual accuracy over an experimental dataset, as parameter re-calibration can signicantly improve the performances of the reduced model on real genetic populations.

Ultimately, the existence of a reduced model will considerably speed up the integration of genetic control into ecophysiological models. Currently, most genetic-improved ecophysiological models make use of Quantitative

Trait loci (QTL) to describe the genetic architecture of specic model parameters. Basically, each parameter has a specic distribution in the population of genotypes and QTL analyses can be performed for each parameter to decipher the architecture of its genetic control (QTL number and eects, linkage).

However, a major drawback of this approach is the diculty in the calibration of the models for a large number of genotypes (due to a large number of parameters along with restricted number of observations) [START_REF] Martre | Chapter 14 -model-assisted phenotyping and ideotype design[END_REF][START_REF] Bertin | Under what circumstances can process-based simulation models link genotype to phenotype for complex traits? case-study of fruit and grain quality traits[END_REF]. Indeed, the statistical power of QTL analyses strongly depends on the size of the population and on the QTL eects i.e. their contribution to the variation of the trait they are associated with [START_REF] Mackay | The genetics of quantitative traits: challenges and prospects[END_REF]. So, in order to be of interest, genetic parameters have to vary among genotypes and be quantiable with relevant accuracy either experimentally or through numerical optimization.

In this perspective, a reduced model with a simpler structure will allow for a better exploration of the parameter space and a more accurate estimation of parameter values. Moreover, the improved calibration time opens the possibility of exploring larger genetic populations so to get more robust QTLs estimation. Finally, it will allow to do simulations over a large number of environmental conditions and/or climatic scenarios.

This is an important step towards dealing with complex Genotype x Environment x Management interactions issues expected in the near future. The development of reliable gene-to-phenotype models will be an important lever to optimize farming in the future climatic conditions.

I(t) = σ f dDW dt + R(t) = (σ f + q g ) dDW dt + q m DW Q (T -20) 10 10 R(t) = q m DW Q (T -20) 10 10 + q g dDW dt DW = DW (t 0 ) + w 1 (1 -e -w 2 t ) + w 3 1+e -w 4 (t-w 5 ) u 1 (I) = 1 σ 1 V 1 λλ suc (t)I(t) λ suc (t) = p 1 t tmax
where t max corresponds to the maturation time The original model by Desnoues et al. [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] was composed by a network 803 of 19 reactions and one input function I(t). The latter described the carbon supply from the mother plant to the fruit and it was estimated as the sum of the carbon used for fruit dry mass (DW ) increase and the carbon lost by respiration (R(t)). Two parameters λ and λ suc described the fraction of the input ow that is converted into the dierent forms of sugars. Fruit respiration was computed following the growth-maintenance paradigm, as described in [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF].

u 2 (I) = 1 σ 2 V 1 (1 -λ)I(t) u 3 (I) = 1 σ 3 V 1 λ 2 (1 -λ suc (t))I(t) Metabolism u 9 (v 2 , x 1 ) = v 2 (t) p 5 +x 1 x 1 (t) u 9 (x 1 ) = v 2 p 5 x 1 (t) = r 1 x 1 (t) u 10 (x 1 ) = v 3 p 21 +x 1 x 1 (t) u 10 (x 1 ) = v 3 p 21 x 1 (t) = r 2 x 1 (t) u 11 (v 4 , x 2 ) = v 4 (t) p 22 +x 2 x 2 (t) u 11 (x 2 ) = v 4 p 22 x 2 = r 3 x 2 (t) u 12 (v 5 , x 2 ) = v 5 (t) p 13 +x 2 x 2 (t) u 12 (x 2 ) = v 5 p 13 x 2 (t) = r 4 x 2 (t) u 13 (x 6 , x 8 , x 9 ) = v 1 (1+
Reaction rates are reported in Table A.3. Enzymatic reactions were generally described using an irreversible Michaelis-Menten kinetics, with experimentallymeasured capacities v i (t). Transport processes between cytosol and vacuole were assumed proportional to the vacuole surface (hypothesis of constant density of transporters) computed from vacuole fresh mass (proxy of the volume) supposing the vacuole as a sphere of surface S(t) = (4π)

1 3 (V 2 ) 2 3
(see [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] for more information). Both active and passive transport mechanisms were considered for fructose and glucose.

Model equations are reported in

dt = u 1 + σ 5 σ 1 u 16 -u 10 -u 4 -µ(t)x 1 dx 2 dt = u 2 -u 11 -u 12 + 1 σ 2 V 1 u 19 -µ(t)x 2 dx 3 dt = u 3 + 1 σ 3 V 1 u 8 + 1 2 σ 1 σ 3 u 9 + 1 2 σ 1 σ 3 u 10 + σ 2 σ 3 u 11 -u 7 -u 14 -µ(t)x 3 dx 3 dt = u 3 + 1 σ 3 V 1 u 8 + 1 2 σ 1 σ 3 u 9 + 1 2 σ 1 σ 3 u 10 + σ 2 σ 3 u 11 -u 14 -µ(t)x 3 dx 4 dt = u 3 + 1 σ 4 V 1 u 6 + 1 2 σ 1 σ 4 u 10 + σ 2 σ 4 u 12 -u 5 -u 15 -µ(t)x 4 dx 4 dt = u 3 + 1 σ 4 V 1 u 6 + 1 2 σ 1 σ 4 u 10 + σ 2 σ 4 u 12 -u 15 -µ(t)x 4 dx 5 dt = 1 2 σ 1 σ 5 u 9 + σ 3 σ 5 u 14 + σ 4 σ 5 u 15 -u 17 -u 16 -1 σ 5 V 1 u 18 -µ(t)x 5 x 5 = 1 p 6 +p 7 +µ(t) ( 1 2 σ 1 σ 5 u 9 + σ 3 σ 5 u 14 + σ 4 σ 5 u 15 -1 σ 5 V 1 u 18 ) dx 6 dt = αu 4 -u 13 -µ(t)x 6 dx 7 dt = -1 σ 7 V 2 u 19 -µ(t)x 7 dx 8 dt = αu 7 + 1 2 σ 6 σ 8 u 13 -1 σ 8 V 2 u 8 -µ(t)x 8 dx 8 dt = 1 2 σ 6 σ 8 u 13 -1 σ 8 V 2 u 8 -µ(t)x 8 dx 9 dt = αu 5 + 1 2 σ 6 σ 9 u 13 -1 σ 9 V 2 u 6 -µ(t)x 9 dx 9 dt = 1 2 σ 6 σ 9 u 13 -1 σ 9 V 2 u 6 -µ(t)x 9 dx 10 dt = σ 5 σ 10 u 17 -µ(t)x 10

Figure 1 :

 1 Figure 1: Graphical representation of the proposed model reduction scheme. Yellow diamonds represent model evaluation steps by means of our 3 criteria: the corrected AIC value, calibration time and expected error over a virtual population. The tested reduction methods are indicated in green. Multivariate sensitivity analysis and three structural simplication methods are independently applied to the original model and evaluated. The validated methods are then combined into an intermediate reduced model whose performances are again submitted to evaluation. Finally, the application of a QSS approximation over the intermediate reduced model is tested to yield the nal reduced model.

Figure 2 :

 2 Figure 2: S 4 is the glucose in the cytosol transported to the vacuole as S 9 via an active (unidirectional transport) and passive (reversible transport).

Figure 3 :

 3 Figure 3: Principal component analysis (PCA) for the whole progeny of 106 genotypes. It represents the projection on the Dim1 and Dim2 of the growth duration and growth rate obtained with growth curves.

S

  : set of the N S simulation times for the genotype k • T (k) M : set of the N M measurement times for the genotype k • X (k) (t j ): N experimental observations for the genotype k, with t j ∈ T (k) M . Note that N = 4 × N M × r, where r is the number of replicates at time t j , for the 4 dierent sugars (sucrose, glucose, fructose and sorbitol). r = 3 for the training set and r = 1 for the validation set.

Figure 4 :

 4 Figure 4: Aggregate Generalized sensitivity indices (aGSI) for the parameters of the model and genotypes (the training set) on four outputs (Sucrose, Sorbitol, Glucose and Fructose) of the sugar model. The main sensitivity indices are in dark bars and interaction ones are in grey bars.

  results shows that the choice of the genotype essentially aects the second principal component, via the denition of the initial conditions of the model (see the supplemental information Fig. B.1). Among the 14 parameters estimated (p 1 , . . . , p 14 ) in the original model, four parameters, namely p 5 , p 10 , p 12 and p 14 , have a negligible eect on the four outputs, independently of the peach genotype. Accordingly, these parameters can be xed to their nominal values i.e. their average value over the ten genotypes, without aecting the quality of predictions. The validity of such a reduction strategy was tested on the ten genotypes of the training set. The dierence in Akaike criterion (∆ AIC C ) between the reduced and the original models was computed for each genotype. Results presented in

Figure 5 :

 5 Figure 5: Dierences in order of magnitude between enzyme anity (K m ) and substrate concentration (x) calculated over the whole dynamics and the training set for each reaction rate u i , i ∈ {5, 7, 9 . . . 15}.

Figure 6 :

 6 Figure 6: Evolution of the active ux (solid lines) and passive transport (dashed lines) for glucose (respectively fructose) and net ux during fruit development (DAB, day after bloom) for the ten genotypes of the training set (dierent colors).

6. 3 .

 3 Strategy 3: Time-scale analysis and QSSA Results from the reduction strategies 1 and 2 were combined into an intermediate reduced model. This model had only 9 parameters to be estimated, linear ows and only one temporal enzymatic capacity, common to all genotypes. Improvement in AIC C with respect to the original model conrmed a strong benece for all ten genotypes (Table

Figure 7 :

 7 Figure 7: Order of magnitude of time scales τ i along fruit development (DAB, days after bloom) for the 10 genotypes of the training set.

Figure 10 :

 10 Figure 10: Evolution of the concentration (mggF W -1 ) of sugars during fruit development (DAB, days after bloom) for ten representative genotypes of the validation set with standard (left) and low fructose (right) phenotypes. Dots represent experimental data and lines are model simulations.

  ), V 2 = max(DAB)/2 + 2, . . . , V 19 = max(DAB)/2 + 19, V 20 = max(DAB))). Results of the principal components and sensitivity principal indices are presented in Fig. B.1. For sucrose, glucose and fructose, the rst two components explained 96% of the total inertia of the simulated sugar dynamics. For sorbitol, the rst three components explained 97%. The rst component was positively correlated with all time-points. Correlation values in Fig. B.1 show that the rst principal component corresponds to the average concentration of sugars (sucrose, glucose, fructose and sorbitol) produced during the whole fruit development. The second principal component was positively correlated with sugar concentration at stage 1 and poorly or slightly negatively correlated with simulated sugar during the second half of fruit development. Thus, the second principal component corresponds to the dierence in sugar initialization values, that strongly depends on the genotype factor. For sorbitol, the third principal component accounts for a much smaller part of inertia, associated with the dierence between the sorbitol produced in the middle of fruit development and the sorbitol produced both very early and late. It was sensitive to the set of genotypes.

Figure B. 2 :

 2 Figure B.2: Generalized sensitivity indices (GSI) for the rst ten sensitive parameters (p i ) and ten genotypes (the training set) on four outputs (Sucrose, Sorbitol, Glucose, and Fructose) of the sugar model. The main sensitivity indices are in dark bars and interaction ones are in grey bars.

Figure C. 3 :

 3 Figure C.3: Principal component analysis (PCA) for the whole progeny of 106 genotypes (grey) and 500, out of 20000, virtual genotypes (black). Represents the projection on the Dim1 and Dim2 of the growth duration and growth rate obtained with curves growth.

Figure D. 4 :

 4 Figure D.4: Evolution of the concentration (mggF W -1 ) of x 5 : Hexose P hosphate during fruit development (DAB, days after bloom) for ten genotypes for the original, intermediate reduced and nal models.

  7 is the vector of time-dependent measured enzymatic activities; p = (p 1 , . . . , p 23 ) is the vector of parameters dening the rate reactions where p 1 , . . . , p 14 have to be estimated and p 15 , . . . , p 23 are xed from literature

	data. f (x(t), I(t), v(t), p) of Eq.(1) describes the change in compounds con-
	centrations. Equations of the reduced and original model are introduced in
	Appendix Appendix A.1.

Table 1 :

 1 Characteristics of enzymatic activities in [1]V max Phenotype eect Temporal eect v 1

		No	No
	v 2	No	Yes
	v 3	Yes	No
	v 4	No	Yes
	v 5	No	Yes
	v 6	Yes	Yes
	v 7	Yes	Yes
	3.2.2. Rate simplication

In the original model, enzymatic reactions were represented by an irreversible Michaelis-Menten (MM) equation:

Table 2

 2

show that such a reduction in the number of parameters is strongly benecial for nine out of the ten genotypes with largely negative ∆ AIC C values, and roughly neutral for one genotype (∆ AIC C ∼ 0). The gain in calibration time, however, is important (25%) and the expected error over the progeny of virtual genotypes is low, demonstrating a good reliability of the proposed simplication. For these reasons, the model with 10 parameters to be estimated was selected.

Table 2 :

 2 ∆ AIC C calculated between reduced and original models for the training set and the gain in calibration time (%) for E43. The Expected error ± standard deviation (Std) between original and dierent reduced models for 20 000 virtual genotypes.

	Simplication method Low sensitive parameters xed Vmax Type eect removed	v3 v6 v7	∆AIC C E1 E33 E43 F111 E22 F106 F146 H191 C216 C227 < ∆AIC C >G Time gain % Calibration -11.5 -6.4 -0.9 -14.04 -13.2 -28.3 -13.5 -14.3 -18.7 -87.7 -20.8 25.8 -1.01 -5.9 -4.15 -4.2 1.1 -2.3 -0.3 -6.1 -6.1 -72.02 -7.9 22.4 -0.1 -4.3 0.06 -3.9 0.7 -5.5 -0.3 -5.4 -6.1 -87.7 -11.3 26.6 -0.7 -36.4 0.2 -6.02 1.4 -5.6 1.9 -5.2 -4.9 -94.5 -14.9 33.9	Expected Error Virtual genotypes 4.9 ± 6.5 0.5 ± 1.3 1.7 ± 1.6 2.9 ± 4.5
	Vmax	v2	-0.1	-3.1	0.06	-3.7	0.7	-5.1	-0.3	-0.3	-6.3	-83.4	-10.1	31.6	0.3 ± 0.7
	Temporal eect removed	v4 v5 v6	-0.8 0.2 -0.3	-8.2 -6.8 -0.4	0.7 0.5 -0.1	-5.6 -4.8 -27.04	-5.03 1.8 1.7	-2.5 -5.8 -5.5	-2.5 2.03 -0.2	-5.1 -2.9 -5.3	-6.1 -6.1 -5.7	-90.3 -91.1 -84.9	-12.3 -11.3 -12.7	19.4 20.3 30.5	2.9 ± 2.5 5.5 ± 5.7 4.1 ± 3.1
		v7	8.6	-25.1	21.1	11.02	19.6	20.01	29.05	12.4	15.6	-97.5	1.5	24.2	6.8 ± 4.5
	Rate simplication Futile cycle removal Intermediate reduced model model Final reduced	Eq.(8) Eq.(9)	-17.2 2.5 0.7 -32.7 -32.5	-53.4 -0.9 -56.7 -18.6 -19.1	8.9 15.6 -5.6 -3.7 -4.3	-35.4 -1.6 -37.1 -24.5 -25.1	-2.9 -0.01 -9.02 -11.8 -12.7	2.7 -2.3 -10.5 -24.04 -1.01-	-14.7 -0.6 -6.7 -16.5 -16.4	-22.9 -1.5 -35.5 -20.3 -20.4	-5.7 -5.9 -12.2 -18.8 -18.8	-71.04 -43.23 -70.7 -43.1 -43.3	-21.1 -3.8 -24.3 -21.4 -18.5	6.7 23.6 24.1 30.5 43.3	18.6 ± 9.7 12.7 ± 14.7 11.5 ± 9.9 22.5 ± 8.4 22.5 ± 8.5
	6.2. Strategy 2: Structural simplication of the model
	Structural simplication methods are another way to reduce the com-
	plexity of dynamic systems by improving the generality of the model and the
	numerical integration of the ordinary dierential equations.

Table A

 A

	.3: Reaction rates of the original and reduced models
	Equations	Original Model	Reduced Model
	Input ows

  ) = r 5 x 6 (t) u 14 (v 6 , x 3 ) = v 6 (t)p 3 +x 3 x 3 (t) u 14 (x 3 ) = v 6 p 3 x 3 (t) = r 6 x 3 (t) u 15 (v 7 , x 4 ) = v 7 (t) p 4 +x 4 x 4 (t) u 15 (v 7 , x 4 ) = v 7 (t) p 4 x 4 (t) u 16 (x 5 ) = p 7 x 5 (t) u 16 (x 5 ) = p 7 x 5 (t) u 17 (x 5 ) = p 6 x 5 (t) u 17 (x 5 ) = p 6 x 5 (t) (S, x 4 , x 9 ) = (x 9 -x 4 )p 10 S(t) u 6 (S, x 4 , x 9 ) = (x 9 -x 4 )p 10 S(t) u 7 (S, x 3 , x 4 ) = (S, x 3 , x 8 ) = (x 8 -x 3 )p 9 S(t) u 8 (S, x 3 , x 8 ) = (x 8 -x 3 )p 9 S(t) u 19 (S, x 2 , x 7 ) = p 14 (x 7 -x 2 )S(t) u 19 (S, x 2 , x 7 ) = p 14 (x 7 -x 2 )S(t)

		x 8 +x 9 p 2 u 13 (x 6 u 18 (R) = R(t) )p 23 +x 6 x 6 (t) u 18 (R) = R(t)
		u 4 (S, x 1 ) = p 8 x 1 (t)S(t)	u 4 (S, x 1 ) = p 8 x 1 (t)S(t)
		u 5 (S, x 3 , x 4 ) =	p 11 p 19 +x 3 +x 4 x 4 (t)S(t)	u 5 = 0
	Transport processes	u 6 p 12 p 20 +x 3 +x 4 x 3 (t)S(t)	u 7 = 0
		u 8

Table A

 A

	.4: System of original and reduced models System of original model	System of reduced model
	dx 1

Table A .

 A 5: Table of original and reduced models parameter description -2.3t + 2e -2 t 2 -8.3e -5 t 3 mg gFW -1 day -1

	pi	Parameter	Corresponding model	Description	Reference	Value	Unit
	p1	λSuc	original and reduced	sucrose proportion hydrolyzed in the apoplasm	Estimated	0-1
	p8	T actif Suc	original and reduced	coecicent of sucrose transport (active import) from cytosol to vacuole	Estimated	0-400	mg gFW -1 day -1
	p10	T passif Glu	reduced	coecient of glucose passive transport between cytosol and vacuole and in the opposite direction	Section. 3.1	112.1559	mg gFW -1 day -1
				original		Estimated	0-150
	p9	T passif F ru	original and reduced	coecient of fructose passive transport between cytosol and vacuole and in the opposite direction	Estimated	0-150	mg gFW -1 day -1
	r1 = v2 p5	Rsusy = Vsusy Ksusy	reduced	coecient of the transfer function between sucrose and (fructose + hexoses phosphate) under action of sucrose synthase (susy) enzyme	Section. 3.2.2	1.8809	day -1
	r2 = v3 p21	Rni = Vni Kni	reduced	coecient of the transfer function between sucrose and (glucose +fructose) under action of neutral invertase (ni) enzyme	Section. 3.2.2	95.5875	day -1
	r3 = v4 p22	Rsdh = Vsdh Ksdh	reduced	coecient of the transfer function between sorbitol and fructose under action of sorbitol dehydrogenase (sdh) enzyme	Section. 3.2.2	7.1592	day -1
	r4 = v5 p13	Rso = Vso Kso	reduced	coecient of the transfer function between sorbitol and glucose under action of sorbitol oxydase (so) enzyme	Estimated	0-10	day -1
	r5 = v1 p23	Rai = Vai Kai	reduced	coecient of the transfer function between sucrose and (glucose +fructose) under action of acid invertase (ai) enzyme	Estimated	0-1	day -1
	p2	KiAI	original	inhibitor constant of acid invertase	Estimated	0-10	mg gFW -1
	r6 = v6 p3	Rfk =	Vfk Kfk	reduced	coecient of the transfer function between fructose and hexoses phosphate under action of fructokinase (fk) enzyme	Estimated	0-5000	day -1
	v7 86.2 p4 V hk(t) reduced hexokinase activity (hk) to transfer glucose to hexoses phosphate Section. 3.2.1 Khk original and reduced hexokinase anity Estimated	1-300	mg gFW -1
	p7	ReSyntSuc	original and reduced	coecient of the transfer function between hexoses phosphate and sucrose	Estimated	0-300	day -1
	p6	OthComp	original and reduced	coecient of the transfer function between hexoses phosphate and other compounds	Estimated	450-1500	day -1
	p14	T passif Sor	reduced	coecient of sorbitol passive transport between cytosol and vacuole	Section. 3.1	4.1305	mg gFW -1 day -1
				original		Estimated	0-150
	σf	P ropCdw	original and reduced	carbon concentration of the mesocarp	[43]	0.44	gC gDW -1
	p17	qg		original and reduced	growth respiration coecient	[43]	0.084	gC gDW -1
	p18	qm		original and reduced	maintenance respiration coecient	[43]	2.76e-4	gC gDW -1 day -1
	p16	Q10		original and reduced	temperature ratio of maintenance respiration	[43]	1.9
	p15	λ		original and reduced	sucrose sap proportion	[1]	0.65
	p11	V mtactif F ru	original	fructose active import (activity)	Estimated	0-150	mg gFW -1 day -1
	p12	V mtactif Glu	original	Glucose active import (activity)	Estimated	0-150	mg gFW -1 day -1
	p19	Kmtactif Glu	original	Glucose active import (anity)	[44]	0.054	mg gFW -1
	p20	Kmtactif F ru

Table D .

 D 6: NRMSE between model simulation and experimental data. Calculated values of the normalized root mean squared error (NRMSE) are presented for each genotype, the four sugars separately.

		Genotype	Phenotype	Year	Sucrose	Sorbitol	Fructose	Glucose	Mean
		E1	Standard	2012	0.09	0.14	0.16	0.26	0.16
		E33	Standard	2012	0.04	0.19	0.27	0.17	0.16
		E43	Standard	2012	0.07	0.11	0.24	0.16	0.14
	Trainig set	F111 C227 E22 F106	Standard Standard Low Low	2012 2011 2012 2012	0.13 0.07 0.11 0.07	0.13 0.28 0.09 0.7	0.21 0.16 0.21 0.14	0.23 0.13 0.18 0.15	0.17 0.16 0.14 0.26
		F146	Low	2012	0.05	0.14	0.02	0.17	0.10
		H191	Low	2012	0.07	0.15	0.18	0.16	0.14
		C216	Low	2011	0.08	0.26	0.25	0.11	0.17
		H163	Standard	2012	0.10	0.35	0.18	0.19	0.20
		F107	Standard	2012	0.11	0.11	0.24	0.25	0.17
		E23	Standard	2012	0.15	0.30	0.10	0.17	0.18
		E17	Standard	2012	0.14	0.37	0.19	0.05	0.18
		E21	Standard	2012	0.13	0.21	0.16	0.24	0.18
		E41	Low	2012	0.08	0.35	0.35	0.24	0.25
		E18	Low	2012	0.13	0.26	0.26	0.09	0.18
		F113	Low	2012	0.12	0.32	0.35	0.20	0.24
		F90	Low	2012	0.06	0.47	0.26	0.13	0.23
		C243	Low	2012	0.18	0.49	0.24	0.09	0.25
		C199	Low	2012	0.22	0.46	0.25	0.19	0.28
		C207	Low	2012	0.19	0.28	0.23	0.24	0.23
	Validation set	E36 E48 F101 F109 F127 F144	Low Low Low Low Low Low	2012 2012 2012 2012 2012 2012	0.09 0.07 0.15 0.07 0.09 0.19	0.13 0.41 0.43 0.31 0.17 0.11	0.13 0.14 0.20 0.27 0.22 0.24	0.16 0.14 0.05 0.24 0.10 0.22	0.12 0.19 0.20 0.22 0.14 0.19
		F141	Low	2012	0.14	0.22	0.36	0.16	0.22
		F86	Low	2012	0.06	0.13	0.30	0.32	0.20
		C232	Standard	2012	0.05	0.30	0.25	0.13	0.18
		E5	Standard	2012	0.22	0.07	0.12	0.24	0.16
		E19	Standard	2012	0.05	0.14	0.24	0.03	0.11
		E20	Standard	2012	0.19	0.17	0.07	0.18	0.11
		E26	Standard	2012	0.05	0.20	0.18	0.38	0.20
		E34	Standard	2012	0.22	0.43	0.20	0.15	0.25
		E35	Standard	2012	0.11	0.27	0.19	0.07	0.16
		E37	Standard	2012	0.24	0.24	0.13	0.26	0.21
		F83	Standard	2012	0.09	0.19	0.06	0.27	0.15
		F115	Standard	2012	0.03	0.20	0.10	0.07	0.10

Acknowledgments

HK was founded by a scholarship of the Lebanese government. We would like to thank V. Signoret for her help in maintaining the inter-specic peach progeny. We are grateful to the IE-EMMAH UMR1114 and IE-GAFL UR1052 teams for taking care of the experimental orchard, and to Dr Olivier Martin for his help in statistics.

Accordingly, x 5 was assumed to be at quasi-steady-state and its equation was replaced by an algebraic function of the slow variables.

Figure 8: Order of magnitude of the predicted sugars concentrations (mg gFW -1 ) in the cytosol (x 1 : Sucrose, x 2 : Sorbitol, x 3 : Fructose, x 4 : Glucose, x 5 : Hexose Phosphate, x 10 : Other compounds) and vacuole (x 6 : Sucrose, x 7 : Fructose, x 8 : Glucose, x 9 : Sorbitol), along fruit development (DAB, days after bloom) for the ten genotypes of the training set.

We compared the intermediate reduced model with its QSS approximation by calculating J i (Eq.( 20)) as explained previously. J i was very low, less than 1%, over the whole dynamics for all variables (Fig. 9). This result was validated also on the virtual genotypes simulated with QQS approximation (see the supplemental information Fig. D.5). In addition the QSS assumption, further increased the performance of the model, leading to a gain in the calibration time of 40% with respect to the original model.

Appendix A.1. Model equations

The original model [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] was written in terms of species carbon quantities C(t). Here, we decided to rewrite the system as a function of species concentration x i (t), for a better readability. The quantity of carbon as a sugar i (C i ) depends on the concentration of i (x i ) according to the following equation:

where σ i is the carbon concentration of sugar i and V j is the volume of the intracellular compartment (cytosol or vacuol) in which species i is located.

The carbon content σ i for the dierent sugar molecules is reported in Table A

Accordingly, for variables 1, . . . , 5, 10,

where x vac i and x cyt i are respectively the variables located in the vacuole (i ∈ [START_REF] Martre | Chapter 14 -model-assisted phenotyping and ideotype design[END_REF][START_REF] Nägele | Mathematical modeling reveals that metabolic feedback regulation of snrk1 and hexokinase is sucient to control sugar homeostasis from energy depletion to full recovery[END_REF]) and cytosol (i ∈ [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF][START_REF] Quilot-Turion | Optimization of allelic combinations controlling parameters of a peach quality model[END_REF]

is the intra-cellular volume ratio. The value of α was estimated by cytological analysis to 0.08 ( see [START_REF] Desnoues | A kinetic model of sugar metabolism in peach fruit reveals a functional hypothesis of a markedly low fructose-to-glucose ratio phenotype[END_REF] for more information). Fruit fresh mass was assumed as a proxy for

For each genotype k, initial conditions x (k) 0 were set using the concentra- tions X (k) (t 0 ) of sucrose, glucose, fructose, sorbitol, and hexose phosphates, measured at the fruit level at the rst stage of development. The conversion between total and intra-cellular metabolite concentrations was performed based on metabolite localization at maturity. Accordingly, 98% of fructose, glucose, sucrose content and 90% of sorbitol content were assumed to be located in the vacuole, whereas the hexose phosphates were restricted to the cytosol. Accordingly, for sucrose, fructose, and glucose:

cytosol :

for sorbitol, cytosol :

and for the hexoses phosphates cytosol : Multivariate sensitivity [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] was used to identify the inuence of each parameter on the dynamic output x(t) during fruit development. Where

x(t) is the sugar concentration (sucrose, glucose, fructose and glucose) and t is the independent time variable for 20 days after bloom (t = (V 1 = Interestingly, the genotype factor is ranked only third to fth, depending on the sugar, meaning that it does not aect parameters' sensitivity as much as

Appendix C. Virtual experiment

In order to evaluate the reliability of the proposed simplications over a larger diversity, a progeny of virtual genotypes was randomly created based on a careful recombination, with noise, of the original 10 dynamics. This included changes in parameters values, initial conditions and input functions.

We used the results from the principal component analysis (PCA) performed on growth rate and growth duration for the whole progeny of 106 genotypes to verify the distribution of virtual genotypes. To this aim, growth rates and durations of the 20 000 virtual genotypes were projected on the PCA plan dened by the previous analysis. As shown in The NRMSE can dened as follows:

where N M is the number of observations, x(t, p(k) ) are the concentrations predicted by the model and X (k) t(t) are the experimental data and i is the sugar index.