M Monneret

Tayakout-Fayolle

J Tallon

B Browning

F Couenne

C Bordes

F Venet

P Nony

F Gueyffier

V Moucadel

G Monneret

M Tayakout-Fayolle
email: melaz.tayakout-fayolle@univ-lyon1.fr

Dynamical Modeling of pro-and anti-inflammatory cytokines in the early stage of septic shock

Keywords: Early sepsis, Dynamic modeling, Parameter estimation, IL18 and IL10 cytokines, Inflammation

Dynamical modeling of pro-and anti-inflammatory

cytokines in the early stage of septic shock

Introduction

Sepsis is defined as a life-threatening organ dysfunction caused by a dysregulated host response to infection [START_REF] Singer | The Third International Consensus Definitions for Sepsis and Septic Shock (Sepsis-3)[END_REF]. Septic shock is the most severe form of sepsis in which hypotension persists despite adequate volume resuscitation thus requiring the use of vasopressors. Sepsis deeply perturbs immune homeostasis by inducing a complex response that varies over time, with the concomitant occurrence of both pro-and anti-inflammatory mechanisms.

Schematically, the opening tremendous systemic inflammatory response (aka cytokine storm) may lead to multiple organ failure while the anti-inflammation process may lead to delayed immunosuppressed status [START_REF] Delano | The immune system's role in sepsis progression, resolution, and long-term outcome[END_REF][START_REF] Hotchkiss | Sepsis-induced immunosuppression from cellular dysfunctions to immunotherapy[END_REF]. However exact chronology of these processes remains unclear. As Sepsis represents a major healthcare problem worldwide (e.g, first cause of mortality in Intensive Care Units (ICU) [START_REF] Sakr | ICON Investigators, Sepsis in Intensive Care Unit Patients: Worldwide Data From the Intensive Care over Nations Audit[END_REF]), the World Health Organization (WHO) recognized sepsis as a global health priority by adopting a resolution to improve the prevention, diagnosis and management of this deadly disease [START_REF] Reinhart | Recognizing Sepsis as a Global Health Priority -A WHO Resolution[END_REF]. Indeed, despite significant advances in the understanding of the pathophysiology of sepsis, to date, no therapeutic intervention targeting host response has specifically been approved.

Several studies have shown that the first hours after septic shock are decisive in the evolution of the disease and therefore for patient's care and outcome [START_REF] Calvano | Lowry and Inflammation and Host Response to Injury Large Scale Collaborative Research Program, A network-based analysis of systemic inflammation in humans[END_REF][START_REF] Cazalis | Early and dynamic changes in gene expression in septic shock patients: a genome-wide approach[END_REF]. At this step, the host response is very complex because a tremendous number of mediators, receptors and cells are involved in the whole body, both activators and inhibitors [START_REF] Yiu | Dynamics of a cytokine storm[END_REF], thus the dynamical modeling approach is an interesting solution to better describe early dynamics in septic shock. In particular, mechanistic modelling based on principles of chemical engineering may help to build the complex models necessary to describe the immune system [START_REF] Rao | On the analysis of complex biological supply chains: From Process Systems Engineering to Quantitative Systems Pharmacology[END_REF][START_REF] Foteinou | Networks, biology and systems engineering: a case study in inflammation[END_REF][START_REF] Aris | Mathematical modeling: a chemical engineer's perspective[END_REF][START_REF] Levenspiel | Chemical reaction engineering[END_REF].

Many dynamical models have appeared in the literature to describe the immune system in the context of sepsis. They are based on ordinary differential equations (ODE) with classical kinetics used in biology such as Michaelis-Menten kinetics, Hill function, inhibitory functions [START_REF] Kumar | The dynamics of acute inflammation[END_REF][START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. Derivation of model and analysis of anti-inflammation[END_REF][START_REF] Day | A reduced mathematical model of the acute inflammatory response II. Capturing scenarios of repeated endotoxin administration[END_REF][START_REF] Scheff | Modeling the influence of circadian rhythms on the acute inflammatory response[END_REF][START_REF] Shi | Mathematical Model of Innate and Adaptive Immunity of Sepsis: A Modeling and Simulation Study of Infectious Disease[END_REF][START_REF] Klett | Modeling the inflammatory response in the hypothalamus ensuing heat stroke: iterative cycle of model calibration, identifiability analysis, experimental design and data collection[END_REF][START_REF] Malek | Dynamics of the HPA axis and inflammatory cytokines: Insights from mathematical modeling[END_REF][START_REF] Wang | The dynamics of NF-κB pathway regulated by circadian clock[END_REF][START_REF] Malkin | A Neutrophil Phenotype Model for Extracorporeal Treatment of Sepsis[END_REF][START_REF] Ma | Defining network topologies that can achieve biochemical adaptation[END_REF][START_REF] Clermont | Using nonlinear model predictive control to find optimal therapeutic strategies to modulate inflammation[END_REF][START_REF] Androulakis | Circadian characteristics of permissive and suppressive effects of cortisol and their role in homeostasis and the acute inflammatory response[END_REF].

A discussion about reduced and large models for immunology and their help with respect to qualitative and quantitative behavior in order to improve diagnosis and treatment is proposed in Vodovotz et al. [START_REF] Vodovotz | Solving Immunology?[END_REF]. Some reduced models are used to analyze the asymptotic behavior of ODE models with respect to key parameters, such as kinetic constants, and show that, depending on parameter values, different immune system scenarios can be represented, as in Kumar et al. [START_REF] Kumar | The dynamics of acute inflammation[END_REF].

Numerous models are used to simulate the immune system for better understanding, [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. Derivation of model and analysis of anti-inflammation[END_REF][START_REF] Scheff | Modeling the influence of circadian rhythms on the acute inflammatory response[END_REF][START_REF] Shi | Mathematical Model of Innate and Adaptive Immunity of Sepsis: A Modeling and Simulation Study of Infectious Disease[END_REF][START_REF] Malek | Dynamics of the HPA axis and inflammatory cytokines: Insights from mathematical modeling[END_REF][START_REF] Ma | Defining network topologies that can achieve biochemical adaptation[END_REF][START_REF] Androulakis | Circadian characteristics of permissive and suppressive effects of cortisol and their role in homeostasis and the acute inflammatory response[END_REF][START_REF] Zhao | Mathematical modeling of septic shock: an innovative tool for assessing therapeutic hypotheses[END_REF][START_REF] Waito | A mathematical model of cutokine dynamics during a cytokine storm[END_REF]. The authors have different modeling approaches, but, currently, the principal limitation is the availability and complexity of experimental data for validation. Some models are used as the foundation for simulations for testing therapy strategies [START_REF] Malkin | A Neutrophil Phenotype Model for Extracorporeal Treatment of Sepsis[END_REF][START_REF] Clermont | Using nonlinear model predictive control to find optimal therapeutic strategies to modulate inflammation[END_REF][START_REF] Zhao | Mathematical modeling of septic shock: an innovative tool for assessing therapeutic hypotheses[END_REF][START_REF] Clermont | In silico design of clinical trials: a method coming of age[END_REF][START_REF] Foteinou | Modeling endotoxininduced systemic inflammation using an indirect response approach[END_REF][START_REF] Bara | Immune therapeutic strategies using optimal controls with L1 and L2 type objectives[END_REF].

Finally, in Yiu et al. [START_REF] Yiu | Dynamics of a cytokine storm[END_REF], based on the 2006 Clinical Trial of TGN1412, the authors propose an eighteenth-order, linear, time-invariant dynamic system to simulate the behavior of nine interacting cytokines based on data obtained from six healthy volunteers that experienced severe inflammatory response during five days. An identification procedure of 90 parameters is performed from cytokine level measurements. The model reveals plausible cause-and-effect relationships among the cytokines showing how each cytokine induces or inhibits others.

For now, the work to model the immune system is only at its beginning. The complexity of the system and the lack of experimental data are the main barriers to advancement. To our knowledge, no article presents a model tested against transcriptomic data.

Microarray-based expression profiling provides an interesting opportunity to gain knowledge on sepsis. This was done in the case of the study on twenty-eight patients in septic shock, called "Sepsis 48h" and conducted in 2009 by the laboratory LCR SEPSIS HCL and bioMerieux [START_REF] Cazalis | Early and dynamic changes in gene expression in septic shock patients: a genome-wide approach[END_REF]. In this study, a blood sample was taken every 6 hours for 48 hours and a sample was taken after 6 days for some patients. For each sample, a complete genomic study was performed.

By virtue of the "Sepsis 48h" databank, it is therefore possible to build an original model that can be validated experimentally against both genetic expression and cell count.

The objective of this study is to propose a pathophysiological dynamical model capable of describing the evolution of inflammation in case of septic shock based on the "Sepsis 48h" data. The chosen cytokines are IL18 and IL10 as well as their associated receptors IL18R and IL10RA. This choice is motivated by the simplicity of the production mechanisms of this pair of cytokines with their associated receptors.

In section 2 we set out the way we selected the training data among the large available dataset.

In section 3 we present the assumptions and the model. In section 4 we describe parameter estimation based on the "Sepsis 48h" data. Finally, in section 5 we discuss the estimated model, the results obtained and we also propose a qualitative validation of the model on the basis of the "Sepsis 48h" data.

Experimental dataset: Sepsis 48h

The "Sepsis 48h" dataset includes clinical characteristics of twenty-eight patients at the onset of septic shock and admitted to two ICUs of a university hospital (see Cazalis et al. [START_REF] Cazalis | Early and dynamic changes in gene expression in septic shock patients: a genome-wide approach[END_REF] for a detailed description of patient characteristics and data treatment). Briefly, a blood sample was taken every 6 hours for 48 hours. The transcribed RNA was extracted from each sample and the gene expressions were measured using microarrays from Affimetrix (GeneChip® Human Genome U133 Plus 2.0). The full data set has not yet been released but data from initial samples and those taken after 24 and 48 hours is available on the GEO DataSets site under accession number GSE57065. So, for each blood sample, data is available is the form of 54675 X-ray fluorescence intensity values. Furthermore, twenty-five control samples have been obtained from healthy volunteers. A statistical comparison of the genetic expressions of samples from healthy and sick patients highlighted 71% of the human genome affected by a septic shock [START_REF] Cazalis | Early and dynamic changes in gene expression in septic shock patients: a genome-wide approach[END_REF].

In addition to this transcriptomic data, a detailed cell count was performed for each sample for the major leukocytes. The main cells of innate immunity were measured: neutrophils, basophils, eosinophils, and monocytes. The main types of lymphocytes were also counted as described in Venet et al. [START_REF] Venet | Early assessment of leukocyte alterations at diagnosis of septic shock[END_REF].

Biological selection of the probesets

The first step of this work was to select the probesets that were related to cytokines and their receptors from the available dataset. First, a global approach was taken, with the selection of all the probesets that were potentially important in describing the variation of the cytokines and their receptors. The probeset references were translated from the microarray output into standard genetic nomenclature found in the literature. For this purpose, the DAVID database (Database for Annotation, Visualization and Integrated Discovery) tool to convert gene identifiers from one type to another was used [START_REF] Reynolds | A reduced mathematical model of the acute inflammatory response: I. Derivation of model and analysis of anti-inflammation[END_REF][START_REF] Yang | A quantitative model of thermal injuryinduced acute inflammation[END_REF].

The genes relating to cytokines and their receptors were then selected with keyword filters: "interleukin"; "TNF"; "IFN"; "TGF" leading to a first reduction of the dataset to 329 probesets of interest. Then, the objective being to describe the inflammation mechanistically and without considering intracellular reactions and intermediates, a deeper analysis of these probesets was made. Data which could not be exploited because it was below the microarray detection limit was removed and, where more than one microarray output referred to the same biological entity, only the most intense signal was selected. Proteins which were unrelated to cytokines and their receptors according to the literature were also removed [START_REF] Genecards® | The Human Gene Database[END_REF]. This allowed selection of 53 probesets related to well-known cytokines and their receptors, all of which are involved in endotoxin tolerance or pathogen recognition.

Here, the strategy is to propose a simplified model based on an immune system of one pair of cytokines with their receptors and the average quantity of leukocytes. It is intended that this be extended in future work. Further consideration of the probesets narrowed down the number of potential cytokines for modelling to 14, each of which has data available for both receptors and cytokine (about 30 probesets) : IL1A, IL1B, IL2, IL3, IL4, IL6, IL10, IL11, IL12A, IL13, IL16, IL18, IL24 and IL27. Of these, IL1, IL4, IL6, IL10, IL18 have been strongly associated with sepsis [START_REF] Cazalis | Early and dynamic changes in gene expression in septic shock patients: a genome-wide approach[END_REF][START_REF] Cai | IL-10 enhances NK proliferation, cytotoxicity and production of IFN-γ when combined with IL-18[END_REF][START_REF] Grobmeyer | Elevation of IL-18 in Human sepsis[END_REF][START_REF] Mommsen | Elevated systemic IL-18 and neopterin levels are associated with posttraumatic complications among patients wit hmultiple injuries: a prospective cohort study[END_REF][START_REF] Fan | Diagnosing sepsis-The role of laboratory medecine[END_REF][START_REF] Lanziotti | Use of biomarkers in pediatric sepsis: literature review[END_REF][START_REF] Dekker | Predictive value of cytokines for developing complications after polytrauma[END_REF][START_REF] Vanden Berge | Paving the way for precision medecin v2.0 in intensive care by profiling necroinflammation in biofluids[END_REF], with 1L1, IL6 and IL18 being pro-inflammatory and IL4 and IL10 anti-inflammatory. From this set of 5 cytokines, IL18/IL10 was identified as the optimal pair for modelling. This cytokine pair is characterized by two cytokines, one with pro-and the other with anti-inflammatory action.

Both IL18 and IL10 cytokines are very often described in the literature dedicated to sepsis [41 -43] and Eidt et al. [START_REF] Eidt | Biochemical and inflammatory aspects in patients with severe sepsis and septic shock: The predictive role of IL-18 in mortality[END_REF] found mortality directly proportional to IL18 plasma levels, which did not occur with other inflammatory mediators whilst Mierzchala-Pasierb et al. [START_REF] Mierzchala-Pasierb | Interleukin-18 serum levels in sepsis: Correlation with disease severity and inflammatory markers[END_REF] found that IL18 can be used to differentiate sepsis and septic shock status better than other biomarkers. IL10 is the anti-inflammatory cytokine. Indeed, the intensities of IL10 and IL18 probesets were strongly detected by the microarrays and significantly increased after septic shock as indicated by t-test (p-value < 0.05, results not shown). The mean X-ray fluorescence signal measurements are given in Figure 1. Moreover, it is well-known that IL10 and IL18 each adsorb onto only one specific receptor [START_REF] Ogata | KEGG: Kyoto Encyclopedia of Genes and Genomes[END_REF]. These receptors were constructed from two proteins (denoted IL10RA and IL10RB for IL10 and IL18R1 and IL18RAP for IL18) that were each described by a single probeset. As IL10RA is specific to IL10 receptor whilst IL10RB is also engaged by other receptors, the variation of IL10RA is followed. IL18R1 and IL18RAP are both specific receptors of the IL18 cytokine [START_REF] Ogata | KEGG: Kyoto Encyclopedia of Genes and Genomes[END_REF]. Here, we considered these two parts as equimolar in the receptor structure.

This allowed the quantity of receptors to be monitored using only the limiting probeset. The selected cytokines and their receptors are presented in Table 1 with their corresponding probeset.

Cytokine/ Cytokine receptor Probeset

IL10 207433_at IL10RA 204912_at IL10RB 209575_at IL18 206295_at IL18R1 206618_at IL18RAP 207072_at
Table 1 : List of the selected probesets related to the selected cytokines and their receptors

In reality, cytokines such as IL18 and IL10 have pro-or anti-inflammatory actions through complex networks involving feedbacks [42, 47 -49]. Here, we are limited to a single cytokine pair. So, prior to the modeling part, the two following assumptions have been made

to describe both the adsorption and production mechanisms of IL10 and IL18 cytokines: i) the cytokines adsorb specifically onto their receptor (IL10RA and IL18R, respectively)

and no interaction occurs with other receptors and ii) the pro-inflammatory cytokine (IL18) stimulates the production of both pro-and antiinflammatory species while IL10 only stimulates the production of anti-inflammatory ones (as reported in Figure 2).

iii) iv) Figure 2 : Schematic representation of adsorption and production mechanisms related to IL18 and IL10 cytokines

Conversion of X-ray fluorescence data into concentration

In order to be able to apply the model to patient data, it was necessary to convert the X-ray fluorescence intensities into concentrations. Indeed, the material balances were made on species in molar quantities.

However, no calibration curve was available. Indeed, the transcribed RNA cannot be rigorously correlated to the amounts of the expressed proteins: many steps subsequent to the formation of the RNA are necessary before obtaining an effector protein. Correlating the X-

Pro-inflammatory cytokine

Anti-inflammatory cytokine

Anti-inflammatory cytokines

Anti-inflammatory cytokines

Pro-inflammatory cytokine ray fluorescence to the concentration is therefore an over-evaluation of the protein concentration in the medium, but is necessary at first.

All the probeset intensities are numerically treated so that they can be compared relative to each other [START_REF] Cazalis | Early and dynamic changes in gene expression in septic shock patients: a genome-wide approach[END_REF]. A calibration curve between X-ray fluorescence intensity and protein concentration can be constructed from measured protein concentrations and the probesets.

During sepsis, protein production is strongly modified: the calibration curve is therefore constructed from control values only. These 22 probesets for healthy individuals make it possible to set maximum ranges of expression values which are correlated with ranges of protein volume concentrations from the literature [START_REF] Kleiner | Cytokine Levels in the Serum of Healthy Subjects[END_REF][START_REF] Sekiyama | Circulating proinflammatory cytokines (IL-1 beta, TNF-alpha, and IL-6) and IL-1 receptor antagonist (IL-1Ra) in fulminant hepatic failure and acute hepatitis[END_REF].

To try and improve accuracy, two other widely measured protein concentrations have been added: S100A8 and S100A9 [START_REF] Cotoi | Plasma S100A8/A9 correlates with blood neutrophil counts, traditional risk factors, and cardiovascular disease in middle-aged healthy individuals[END_REF]. These S100 alarmin biomarkers are secreted by leukocytes and are involved in various inflammatory diseases. This data was used to create the calibration curve shown in Figure 3, relating the protein concentration in mol.m -3 to the X-ray fluorescence intensity. Linear regression was used to obtain the following relation:

Protein concentration = 1.27x10 -9 X-ray fluorescence intensity So, finally, the experimental training data, which will be used for parameter estimation, is the concentrations of IL10, IL10R, IL18, IL18R and the white blood cell count (leukocytes).

Modelling

The model is based on chemical reaction engineering principals with the following assumptions:

 The blood system is an open system (mass balances based on fluxes), described with two phases: a fluid phase (blood), and leucocyte cells. It has interactions with surroundings through source terms of cells coming from bone marrow and source terms of cytokines and receptors coming from organs. The blood system is represented in Figure 3.  The linear driving force model is used to represent the transfer of cytokines from the cell to the blood and the transfer of cytokines from the blood to receptors [START_REF] Bird | Transport Phenomena[END_REF]. This model relates the average adsorbate concentration inside the cell directly with the concentration in the fluid phase.

 The cytokines in the blood are adsorbed on receptors. This adsorption is already described in some existing pharmacokinetic models [START_REF] Androulakis | Circadian characteristics of permissive and suppressive effects of cortisol and their role in homeostasis and the acute inflammatory response[END_REF]. Langmuir equilibrium is assumed for the adsorption [START_REF] Langmuir | The adsorption of gases on plane surfaces of glass, mica, and platinum[END_REF]. So the expressions of the adsorbed concentrations onto the surface of cells are given below (see Nomenclature section for definition of variables):

𝐴 𝑎𝑑𝑠 = 𝑅 𝐴 𝐾 𝐴 𝐴 𝑖𝑛𝑡 1 + 𝐾 𝐴 𝐴 𝑖𝑛𝑡 , 𝐵 𝑎𝑑𝑠 = 𝑅 𝐵 𝐾 𝐵 𝐵 𝑖𝑛𝑡 1 + 𝐾 𝐵 𝐵 𝑖𝑛𝑡
The inflammation comes from the variation of the cytokines, which can stimulate or inhibit leucocyte production.

Cytokine material balances in the fluid, at the cell/fluid interface and in the cells

In chemical reaction engineering it is usual to construct balances over a defined volume based on conservation of mass. This can be for any individual component or the sum of all the species present. The mass balance takes into account the consumption, production and accumulation of the species under consideration as well as mass flows into and out of the defined volume [START_REF] Levenspiel | Chemical Reaction Engineering[END_REF].

The material balances in mol.min -1 for A (pro-inflammatory cytokine IL18) and B (antiinflammatory cytokine IL10) in the fluid, at the interface of the cells and in the cells, are given in Eq.s ( 1) -( 6). Eq. ( 1) is for the pro-inflammatory cytokines in the blood volume.

They arrive from the cell interior and a source, such as an inflamed organ or mucus, and are transferred to receptors on the cell membrane. There is also a term for cytokine consumption because they have a fixed lifetime. Eq.s (2) and (3) respectively, are the balances on the cytokine quantities at the cell membrane and inside the cell. The amount of cytokine interior to the cell depends on the mass transfer rate from the cell to the blood volume, the consumption rate due to the fixed lifetime and the production of cytokine due to adsorption onto the cell surface. Eq.s (4) to ( 6) represent the anti-inflammatory cytokine behaviour. It is identical to that of the pro-inflammatory cytokine except adsorbed quantities of both pro-and anti-inflammatory cytokines are used to determine the cytokine quantity inside the cell. A and B at the cell surface to be managed. For 𝒌 𝟏 and 𝒌 𝟐 equal to 1, A and B receptor concentrations are maintained constant with respect to cell number. For 𝒌 𝟏 and 𝒌 𝟐 greater than 1, A and B receptor densities increase, while for 𝒌 𝟏 and 𝒌 𝟐 lower than 1, they decrease.

In the particular case where 𝒌 𝟏 and 𝒌 𝟐 are equal to 0, the quantity of A and B receptors remains constant.

𝑑 𝑉

Cell number balances

Eq. ( 9) is the balance on the number of cells in the blood volume with the inflammation function given in Eq. [START_REF] Foteinou | Networks, biology and systems engineering: a case study in inflammation[END_REF]. This function combines the quantities of pro-and antiinflammatory cytokines to give a numerical representation of the overall amount of inflammation in the body. This function is chosen with the variable parameter, α, adjusted such that f=0 for healthy volunteers, with the average IL18 and IL10 concentrations taken as

Results and discussion

The estimation procedure was initiated from simulation data reported in the previous work [START_REF] Tallon | Immune system modelling in case of a septic shock[END_REF]. For the pair of cytokines under consideration, the model has 21 parameters for 45 experimental data points per patient (5 components multiplied by 9 sample times) and there are 19 patients. It is thus important to reduce the number of parameters to avoid over parameterization and to this end ten of the parameters (see ) had little impact. Each of these was therefore fixed at an estimated average value. The kinetic constant for cell natural death was also found not to be very sensitive and was fixed to be the same as a typical cell lifetime, around 24h.

Finally, high values were chosen for the mass transfer coefficients so that mass transfer would not be a limiting factor. This seems to be a good assumption because we know the cytokines are transferred rapidly around the whole system. Table 2 : Fixed parameters and their values 11 parameters were estimated from 45 data points using the MATLAB non-linear least squares solver function, lsqnonlin, with the trust-region-reflective algorithm [START_REF]MATLAB version 9.3.0[END_REF]. The algorithm minimizes the objective function based on an input vector of differences between the measured and calculated data.

Parameter significance levels and confidence limits of the vector parameter b were determined from the standard calculation method assuming that errors in the data are normally distributed and bearing in mind the nonlinearity of this model. The lsqnonlin Jacobian output matrix, J, was used to estimate the standard error 𝒔(𝒃 𝒊 ) for parameter 𝒃 𝒊 from Eqs. ( 11) to [START_REF] Levenspiel | Chemical reaction engineering[END_REF] with 𝚺, n and p the sum of the squared errors, the numbers of data points and the number of parameters to estimate respectively. 𝑯 𝒊𝒊 𝑻 is the i th diagonal element of 𝑯 𝑻 . Hence the statistical significance of each parameter in Eq. ( 13) was used to calculate the confidence limits at 𝜶 𝐭 % using the t-distribution with variable t.

𝐻 = 𝐽 𝑇 𝐽 (11) 𝑠(𝑏 𝑖 ) = √ Σ. 𝐻 𝑖𝑖 𝑇 (𝑛 -𝑝) (12)
𝑏 𝑖 ± 𝑠(𝑏 𝑖 )𝑡(𝑛 -𝑝; 𝛼 t /2) (13)

Strategy for parameter estimation

Although 10 parameters have been fixed due to their low variability, 11 parameters (8 physical parameters and 3 source terms) remained to be estimated for a total of 45 data points.

First the estimation of all these parameters was performed for each patient. From a global analysis of parameter sensitivity, the three least sensitive parameters were then fixed and a second round of estimation was performed for all patients. Then, three more parameters were fixed and a third, final, estimation was carried out. In this section, we present the results of the different steps of this strategy.

First round of estimation: 11 parameters for each patient

The estimated parameters are listed in Table 3 and Table 4. The distributions by number of patients for all the estimated parameters listed in Table 3 are given in Figure 4. For each parameter, the estimated values for all the patients were collated into 5 sets to smooth the results and highlight the overall trend. Figure 5 gives the estimated values of the cytokine source terms, 𝑆 𝐴 and 𝑆 𝐵 , plotted against the estimated cell source term, 𝑆 𝑐 . This figure shows that the production of cytokines is more or less independent of the cell source term, except for the pro-inflammatory cytokine source term 𝑆 𝐴 which seems to decrease as the cell source term increases. The bone marrow of women appears to produce fewer cells than that of men. Moreover, Figure 5 shows that the relative rates of cytokine and cell production do not differ between patients who survive and those who do not. The mean parameter values (corresponding to the first moment of the distribution) are listed in Table 5 with parameter values and 95% confidence intervals for four patients (2 men and 2 women including 2 patients who died). The source terms of cytokines and cells are presented in Table 6. (1 woman and 3 men) with their 95% confidence interval (CI).

𝐾 𝐴 𝐾 𝐵 1 k A 1 k B 1 k BA k 1 k 2 1 k c
Confidence intervals were calculated for the nonlinear parameters. A wide confidence interval suggests that there is insufficient identifiability structure in the model to determine the parameters from the available measurements. The existence of superfluous parameters in the model may lead to a "rank deficient" condition of the Jacobian matrix (when gradient based methods are used for solution) and/or inflated confidence intervals.

From the results presented in Table 5, it can be seen that the confidence intervals of parameters k B , k BA , k c were very wide. So, they were considered non-sensitive and were set at the mean value obtained from all patients and given in Table 5. The most sensitive parameter for all patients, without exception, is the source term 𝑺 𝑪 corresponding to the production of the cells by the bone marrow.

4.1.2. Second round of estimation: 8 parameters for each patient

The same methodology as in section 4.1.1 was applied for this new estimation. The distributions of parameters 𝑲 𝑨 , 𝑲 𝑩 , k A , k 1 , k 2 by number of patients are presented in Figure 6. The estimated cytokine source terms, 𝑆 𝐴 and 𝑆 𝐵 , considered relative to the cell source term, 𝑆 𝑐 , are not given because the results were the same as those shown in Figure 5. Table 7 gives the Since the confidence interval of parameter K B is wide for most of the patients, it was set to the corresponding total mean value given in Table 7. From Figure 6 and Table 7, it can be seen that the mean value of parameter k 2 was close to 1 for most patients. The concentration of the associated receptor (R B ) was constant, meaning that the receptor concentration is linked to the formation and death of cells. The variation of the cell number was not accompanied by a proportional change in the number of receptors. So, this parameter was fixed at one for the last estimation. The source term related to the anti-inflammatory components was less sensitive and remained relatively constant for all patients (see Figure 5). This latter was thus also fixed to 0.0742 mol/min, i.e. the average value calculated across all patients. Table 8 gives the mean parameter values and source terms and also the parameter values and source terms for the same four patients as previously with their confidence intervals. (mean values, n=19) and individually for four patients (1 woman and 3 men) with their 95% confidence interval (CI).

Figure 8 shows that the parameters of all patients who died are to the left of the distribution, especially for the adsorption of the pro-inflammatory cytokine, IL18, onto its receptor, meaning a very slow adsorption. This slow adsorption could cause a slow cell response to the inflammatory action of the system and therefore lead to bad regulation of this aspect.

Figure 9 shows the same trends as Figure 5 with the pro-inflammatory cytokine source term decreasing as the cell source term increases for females and no differentiation between surviving and deceased patients.

Model outputs

The parity plots in Figure 10 However, there is no common trend between the four cases shown. The profiles of the curves relating to the women show increasing quantities of IL18 which would be expected at the start of sepsis [START_REF] Yiu | Dynamics of a cytokine storm[END_REF], whereas the curves for the men have much less of a gradient. One drawback with the data is that the first measurement is based on the moment the patient arrived in the hospital and not the actual onset of sepsis. So the initial condition for the model is at an unknown time during the sepsis response. Another is that although the treatment is standardized the cause of sepsis is not controlled.

Living woman Dead Man (1)

Living man Dead man (2) Figure 12 shows the results for the pro-inflammatory cytokine receptor. In all cases, except that of the surviving man, the model represents the data well, suggesting that the method of determining receptor concentration from the quantity of cytokine at the cell membrane is a useful one. Also, the number of data points is quite limited and the timing of the data collection, in the first 48h of sepsis, is during a period of strong inflammation. This could partially explain why the model was more sensitive to parameters relating to the proinflammatory cytokine, IL18, and why the calculated fluorescence for IL18 and its receptor fit the measured data better than those for IL10. Furthermore, in their work to model the response of 9 cytokines to TGN1412 infusion, Yiu et al. [START_REF] Yiu | Dynamics of a cytokine storm[END_REF] demonstrated that IL10 concentrations have a small but rapid response to IFN-γ stimulus and, in their model of cytokine dynamics during a cytokine storm in mice, Waito et al. [START_REF] Waito | A mathematical model of cutokine dynamics during a cytokine storm[END_REF] showed that IL10 concentrations depend on at least six pro-inflammatory cytokines. So, as more data becomes available and more cytokines are included in the reaction network, calculation of the antiinflammatory response should improve.

Living woman Dead Man (1)

Living man Dead man (2) Figure 15 shows that the model gives a good prediction of the variation of the number of leukocyte cells, particularly for the women in this case. This depends on the parameters: k d , which fixes the cell lifetime, S c , which determines the quantity of leukocytes entering the blood volume and , the inflammation function. Figure 16 shows the calculated and measured inflammation level function 𝒇 versus time for the four patients. This inflammation level doesn't provide any general information about the inflammatory state of the patient but only describes the amounts of A (pro-inflammatory cytokine IL18) and B (anti-inflammatory cytokine IL10) cytokine pair. The experimental data is dispersed making its interpretation difficult. The simulation showed increasing inflammation in all cases except for the surviving man and a greater increase in the inflammatory level of the women relative to the men. In fact, the curves for the women are highly reminiscent of those used by Hotchkiss et al. [START_REF] Hotchkiss | Sepsis-induced immunosuppression from cellular dysfunctions to immunotherapy[END_REF] to describe the host inflammatory response in their competing theories of the host immune response in sepsis, with a rapid initial increase in inflammation over the first two days which then plateaus before decreasing.

Living woman Dead Man (1)

Living man Dead man (2)

Conclusion

A dynamical model describing the physico-chemical phenomena involved in homeostasis and sepsis, has been proposed in Tallon et al. [START_REF] Tallon | Immune system modelling in case of a septic shock[END_REF]. This dynamical model of the blood system considers interactions between cells and cytokines, with its surroundings (organs, bone marrow) and introduces the inflammatory level impact on homeostasis via the function f.

Data for the cytokine pair IL18 and IL10 from the sepsis 48H data base has been chosen to validate this model. Since the model contains 21 parameters, 10 parameters (kinetic constants and source terms) have been fixed and the others have been estimated in three stages. All estimated kinetic constant values are distributed with respect to the number of patients.

Finally, the more sensitive parameters are those related to pro-inflammatory cytokines and cells. By the parity diagram, the model shows a good agreement between the experimental and simulated data. At the beginning of sepsis, the simulated patient response seems to be driven by the inflammatory actors. The estimated adsorption coefficient of the proinflammatory cytokine, IL18, on its associated receptor, IL18R, is consistently below the mean value of its distribution and, is particularly low in the non-surviving patients, hinting that this could be a key aspect in the dysregulation of the immune system in sepsis and would be interesting to study with a larger data set.

Nomenclature

Figure 1 :

 1 Figure 1: Mean X-ray fluorescence signal measured for IL18 and IL10 cytokines at 0, 24h and

Figure 3 :

 3 Figure 3: Calibration curve for protein concentration from Average X-ray fluorescence values



  The blood volume is assumed to be constant. The cells are uniformly dispersed in the fluid phase and represent average leucocyte cells. The receptors are described as being only on the surface of the cells. They are said to be uniformly distributed. The adsorbed anti-inflammatory IL10 cytokines (B) on their receptors IL10RA (𝑅 𝐵 ) directly trigger the chemical production of IL10 (B) and IL10RA (𝑅 𝐵 ). The adsorbed pro-inflammatory IL18 cytokines (A) on their receptors IL18R ( 𝑅 𝐴 ) trigger the chemical production of all cytokines and receptors. The production kinetics of both cytokines and receptors are assumed to be of order 1 with respect to the adsorbed cytokines.

2 .

 2 𝑐 (𝑘 𝑙𝑎 𝑖𝑛𝑡𝐴 (𝐴 𝑖𝑛𝑡 -𝐴) + 𝑘 𝑙𝑎 𝑐𝐴 (𝐴 𝑐 -𝐴)) -𝑉 𝑙 𝑘 𝑑 𝐴 𝐴 + 𝑆 𝐴 (1) Accumulation of A in blood (mol.min -1 ) Transfer of A from blood to cell surface Transfer of A from inside cells to blood Consumption of A in blood A into blood from source 𝑑𝑉 𝑐 𝐴 𝑎𝑑𝑠 𝑑𝑡 = -𝑉 𝑐 𝑘 𝑙𝑎 𝑖𝑛𝑡𝐴 (𝐴 𝑖𝑛𝑡 -𝐴) (2) 𝑑𝑉 𝑐 𝐴 𝑐 𝑑𝑡 = 𝑉 𝑐 𝑘 𝐴 𝐴 𝑎𝑑𝑠 -𝑉 𝑐 𝑘 𝑙𝑎 𝑐𝐴 (𝐴 𝑐 -𝐴) -𝑉 𝑐 𝑘 𝑑 𝐴 𝐴 𝑐 (3) 𝑑𝑉 𝑙 𝐵 𝑑𝑡 = 𝑉 𝑐 (𝑘 𝑙𝑎 𝑖𝑛𝑡𝐵 (𝐵 𝑖𝑛𝑡 -𝐵) + 𝑘 𝑙𝑎 𝑐𝐵 (𝐵 𝑐 -𝐵)) -𝑉 𝑙 𝑘 𝑑 𝐵 𝐵 + 𝑆 𝐵 (4) 𝑑𝑉 𝑐 𝐵 𝑎𝑑𝑠 𝑑𝑡 = -𝑉 𝑐 𝑘 𝑙𝑎 𝑖𝑛𝑡𝐵 (𝐵 𝑖𝑛𝑡 -𝐵) (5) 𝑑 𝑉 𝑐 𝐵 𝑐 𝑑𝑡 = 𝑉 𝑐 𝑘 𝐵 𝐵 𝑎𝑑𝑠 + 𝑉 𝑐 𝑘 𝐵𝐴 𝐴 𝑎𝑑𝑠 -𝑉 𝑐 𝑘 𝑙𝑎 𝑐𝐵 (𝐵 𝑐 -𝐵) -𝑉 𝑐 𝑘 𝑑 𝐵 𝐵 𝑐 (6) Accumulation of A at cell surface (mol. min -1 ) Transfer of A from blood to cell surface Transfer of A from inside cells to blood Production of A in cell Accumulation of A in cell (mol. min -1 ) Consumption of A in cells Accumulation of B in blood (mol. min -1 ) Receptor material balances at the cell/fluid interface Eq.s (7) and (8) give the mass balances on the cytokine receptors which are found on the outside surface of the cell membrane. The terms 𝒌 𝟏 𝑽 𝟏𝒄 𝑽 𝒍 𝒅 𝑵 𝒄 𝒅𝒕 𝑹 𝑨 and 𝒌 𝟐 𝑽 𝟏𝒄 𝑽 𝒍 𝒅 𝑵 𝒄 𝒅𝒕 𝑹 𝑩 respectively allow the receptor density of

3 .

 3 87 and 4.09 nmol.m -3 respectively. 𝑑 𝑉 𝑙 𝑁 𝑐 𝑑𝑡 = 𝑉 𝑙 𝑓 𝑘 𝑐 𝑁 𝑐 -𝑉 𝑙 𝑘 𝑑 𝑐 𝑁 𝑐 + 𝑆 𝑐 (9) 𝑓(𝑡) = 𝐴 -𝛼 * 𝐵 (10) This model has the ability to evolve to an alternative homeostatic equilibrium in the case of septic shock as shown in Tallon et al. [56].

Figure 4 : 23 Figure 5 :

 4235 Figure 4: Representation of the blood system

Figure 6 :

 6 Figure 6: Cytokine source terms, 𝑆 𝐴 and 𝑆 𝐵 , relative to cell source term, 𝑆 𝑐 . Red dots indicate

  and source terms and also the parameter values and source terms for the same four patients as previously with their confidence intervals.

Table 7 : 4 . 1 . 3 . 29 Figure 7 :

 7413297 Figure 7: Parameter value distributions by number of patients for 𝑲 𝑨 , 𝑲 𝑩 , k A , k 1 , k 2 .

Figure 8 : 32 Figure 9 :

 8329 Figure 8: Parameter value distributions by number of patients for 𝑲 𝑨 , k A , k 1 and the intervals corresponding to all dead patients. Overall, the parameter values are of the correct order and comparable with other values found

Figure 10 :Figure 11 :

 1011 Figure 10: Parity plots for calculated N c , A, B, R A , R B fluorescence against measured data (each color corresponds to one patient). Pearsons correlation coefficient values are 0.906 for N c , 0.715 for A, 0.818 for B, 0.855 for R A and 0.874 for R B .

Figure 11 compares

 11 Figure 11 compares the measured and calculated values of fluorescence for IL18, the proinflammatory cytokine. In each case, the calculated values follow the trend of the data.

Figure 12 :

 12 Figure 12: Calculated and measured recepetor of pro-inflammatory cytokine R A fluoresence versus time, markers indicate measured data, simulation results are shown as lines.

Figure 13 Figure 13 :Figure 14 :

 131314 Figure 13: Calculated and measured anti-inflammatory cytokine B fluoresence versus time, markers indicate measured data, simulation results are shown as linesthe anti-inflammatory cytokine.As mentioned earlier, the model was not very sensitive to the parameters for the anti-inflammatory cytokine and here we observe that, despite apparent trends in the data, the model did not pick up this dynamic but rapidly adjusted to a value close to the mean.

Figure 14

 14 Figure 14 compares the measured and calculated values of the anti-inflammatory cytokine receptor fluorescence and, here, the receptor of anti-inflammatory cytokine R B does not vary much experimentally and the calculated values were constant. This low variation observed experimentally is almost definitely the reason behind the low-sensitivity of the associated parameters (K B, 𝐤 𝐁 , 𝐤 𝐁𝐀 ). Also, the number of data points is quite limited and the timing of the

Figure 15 :

 15 Figure 15: Calculated and measured cell number N c versus time, markers indicate measured data, simulation results are shown as lines.

Figure 16 :

 16 Figure 16: calculated and measured inflammation level function 𝒇 versus time, markers indicate measured data, simulation results are shown as lines.

3 )

 3 Concentration in the blood of cyt. A, B resp. 𝐴 𝑐 , 𝐵 𝑐 (𝑚𝑜𝑙 𝑚 -3 ) Concentration in the cell of cyt. A, B resp. . A, B resp. 𝑓 (𝑚𝑜𝑙 𝑚 -3 ) the inflammation level function 𝑘 𝐴 , 𝑘 𝐵 (𝑚𝑖𝑛 -1 ) Production kinetic constant of cyt. A, B in the cell from adsorbed cyt. A, B resp.

  𝑐 𝑅 𝐴 𝑑𝑡 = 𝑉 𝑐 𝑘 𝑅𝐴 𝐴 𝑎𝑑𝑠 + 𝑘 1 𝑉 1𝑐 𝑉 𝑙 𝑑𝑉 𝑐 𝑅 𝐵 𝑑𝑡 = 𝑉 𝑐 𝑘 𝑅𝐵 𝐴 𝐴 𝑎𝑑𝑠 + 𝑉 𝑐 𝑘 𝑅𝐵 𝐵 𝑎𝑑𝑠 + 𝑘 2 𝑉 1𝑐 𝑉 𝑙

			𝑑 𝑁 𝑐 𝑑𝑡	𝑅 𝐴			(7)
	Accumulation	Production of	A receptor
	(mol. min -1 ) receptors of A	A receptors	density term	𝑑 𝑁 𝑐 𝑑𝑡	𝑅 𝐵	(8)
	Accumulation	Production of	Production of	B receptor
	of B receptors (mol. min -1 )	due to B receptors	due to B receptors	density term
		adsorbed A	adsorbed B

Table 2

 2

	Accumulation	Production	Consumption	Cells into
	of cells	of cells	of cells	blood from
	(number of			source
	leukocytes.
	min -1 )

) have been fixed at reasonable values. These were chosen by sensitivity analysis of the model, which revealed

Table 3 :

 3 List of the estimated parameters

	Source terms for:

Table 4 :

 4 List of the estimated source terms

Table 5 :

 5 Mean parameter values and calculated parameter values for four Patients (1 woman and 3 men) with their 95% confidence intervals (CI).

	Source terms	𝑺 𝑨	𝑺 𝑩	𝑺 𝑪
		(mol.min -1 )	(mol.min -1 )	(min -1 )
	Living woman	0.14	0.10	0.0008
	CI for Living woman	±0.07	±0.04	±0.0002
	Dead man (1)	0.09	0.07	0.0046
	CI for Dead man (1)	±0.40	±0.30	±1 10 -5
	Living man	0.05	0.07	0.0010
	CI *for Living man	±0.11	±0.12	±0.0003
	Dead man (2)	0.10	0.10	0.0014
	CI for Dead man (2)	±0.03	±0.22	±0.0001

Table 6 :

 6 Mean source term values and source term values calculated for four Patients

Table 8 :

 8 Parameter and source term values determined from the whole set of patients

	k 1	S A	S A
	(-)	(mol.min -1 )	(mol.min -1 )


Tzonka Mineva
email: tzonka.mineva@enscm.fr

] Eddy Dib

Arnold Gaje

Hugo Petitjean

Jean-Louis Bantignies

Bruno Alonso
email: bruno.alonso@enscm.fr

Zeolite structure direction: identification, strength and involvement of weak CH…O hydrogen bonds

Keywords: Hydrogen bonds, Zeolites, Structure directing agent, IR spectroscopy, Density functional calculations, Ab-initio molecular dynamics

We demonstrate that weak CH…O hydrogen bonds (HBs) are important host-guest interactions in zeolite assemblies involving structure directing organocations. This type of HB are identified between alkyl groups of the organic structure directing agent (OSDA) and the silica framework in assynthesized silicalite-1 of complex topology (MFI) using a combination of experimental and theoretical data obtained at low and room temperatures. The 28 weak CH…O HB, evidenced along dynamics simulation at room temperature, represent 30% of the energy of the Coulomb electrostatic interaction between OSDA and the zeolite framework. The strongest and most stable HB found here connects the OSDA to the [4 1 5 2 6 2 ] cage containing F atoms and should contribute to preserve zeolite topology during crystal growth. An inspection of other as-synthesized zeolites of very different framework topology indicates that the directional CH...O HBs have to be considered when discussing zeolite structure directing phenomena.

Since the first systematic investigation, [1] weak hydrogen bonds (HBs) between CH donor groups and O acceptor atoms have been evidenced and described for several organic and biological molecules. [2] Today these HBs have a recognized role in chemistry, notably in the understanding of conformational states [3] and in crystal engineering. [2c,4] In the case of the important class of oxide-based materials that are zeolites, the implication of these weak HBs in crystal formation is scarcely investigated despite their possible structure directing and stabilizing effects. For instance tetraalkylammonium cations are widely used as organic structure directing agents (OSDAs) in the hydrothermal synthesis of zeolites, [5] and are also known for a long time to form CH…X HBs with proton acceptor anions. [6] Besides, Behrens et al. showed the presence of HBs between the aromatic CH groups of the organometallic cobaltocenium cation [Co(C 5 H 5 ) 2 ] + and the O atoms of silica in the particular case of a clathrasil framework (nonasil, NON). [7] Some years later, Hong et al. deducted the existence of intramolecular CH…O HBs for choline cations inside UZM-22 zeolites (MEI framework), [8] not involving directly the zeolite framework. Recently, we evidenced the possible stabilizing role of intermolecular CH…O HBs between tetramethylammonium cations and silica inside another clathrasil framework (octadecasil, AST) thanks to a combination of experimental methods (NMR, XRD) and theoretical calculations. [9] Weak CH…O HBs are not easily identified as their individual signature can be very subtle, and often blurred in the ensemble of intermolecular interactions. In the case of OSDA-zeolite framework interactions, the difficulty is enhanced by the static and thermal disorders at the crystallographic positions. Moreover, the large number of unequivalent atoms in zeolites with complex topologies [10] still hinders the identification of CH…O H-bonding. Despite these difficulties, considering the CH…O HBs when discussing OSDA-zeolite interactions would be beneficial for the design of the final zeolite structure and properties. In this letter, we address for the first time the occurrence and characteristics of these HBs in the topologically complex MFI zeolites of high scientific and industrial interest, [5,11] by inspecting geometrical and spectroscopic (IR) data, and evaluating their strength with respect to the total OSDA-zeolite interactions, including these at finite temperature through density-functionaltheory (DFT) methods and Born-Oppenheimer molecular dynamics (BOMD) simulation.

A first evidence of the presence of CH…O bonds comes through the analysis of crystal structures of as-synthesized zeolites with MFI topology containing tetrapropylammonium (TPA). In these X-ray structures, a great number of C…O distances d CO are below 4.0 Å (Supporting Information, Table S1) suggesting the formation of weak HBs between H donors and O acceptors. [2c] There is also a significant number of cases where d CO < 3.6 Å more certainly ascribed to HB interactions, [2b,12] notably for the TPA-F-silicalite-1 structure. [13] This zeolite synthesized in fluoride medium presents well ordered structures with a low level of chemical defects. It is chosen here to explore the existence and role of CH...O HBs.

A second experimental indication of CH…O HBs is obtained from IR spectroscopy. The variations of the CH stretching vibrations (CH) of the OSDA in TPA-F-silicalite-1 at 275 K and 180 K, above the orthorhombic-monoclinic phase transition, [14] are presented in Figure 1. They exhibit the typical (CH) complex pattern of this material. [15] From IR measurements as a function of the temperature (Supporting Information, Figure S3) we observe two significant gradual downshifts for (CH) localized at ~2850 ( 275-180K) = 4cm -1 ) and ~2917 cm -1 ( 275-180K) = 2cm -1 ) assigned respectively to symmetric and antisymmetric modes of CH 2 . [START_REF] Socrates | Infrared and Raman characteristic group frequencies: tables and charts[END_REF] This non classical pseudo-harmonic behavior is commonly reported when the internal stretching mode are affected by H-bonding interactions [START_REF] Bantignies | [END_REF] meaning that at least one type of CH x groups is involved in HBs. Furthermore, such thermal variations are observed for (CH) in cellulose crystals, [18] and have been recently assigned to CH…O HBs. [19] S1). For the two space groups considered -Pnma the most often determined for silicalite-1 containing TPA, and Pn2 1 a that better accommodates the symmetry of TPA and zeolite frameworkwe found a great number of H…O distances d HO below 3.0 Å, and a significant number (9 to 12) for which d HO < 2.6 Å with CHO angles between 100 and 180°. These geometrical parameters are indeed consistent with the presence of weak HBs. [2b,20] To circumvolve possible artifacts generated by DFT optimizations at T = 0 K, we carried out dynamics simulation (for 1.1 ps) at T = 300 K of a cluster model centered on the TPA position. The spatial organization of the silica shell and the motions experienced by TPA appear here highly compatible with the experimental crystallographic and 13 C and 14 N NMR data (Supporting Information S6), and importantly the TPA dynamics can be appropriately used to investigate CH…O interactions. In the resulting time-averaged cluster (Figure 2A), 30 d HO values are falling below 3.0 Å and involve 19 H (over 28 H in TPA) and 23 O (Supporting Information, Tables S1,S7a). There is no specific cycle or cage location for the O atoms that can form one, two or three HBs although the shortest (and strongest, vide infra) HB is found for O in a 4 MR in the [4 1 5 2 6 2 ] cage containing one F atom. The 3.0 Å per CHO angle range after cone correction, [21] show a distinct preference for the highest angle values (170-180°). Besides, CHO angles in the time-averaged cluster tend to decrease when d HO increase, while remaining above 90° (Supporting Information, Figure S7a). Such behaviors reflect the directionality of HBs (stronger when the angle is closer to 180°). [20] They confirm indeed the presence of CH…O HBs allowing to discard van der Waals interactions for these CHO moieties. [12] The dynamics results are even more informative. The d HO profiles as a function of the computed time allow to distinguish between: i) stable HBs (d HO < 3 Å all along the dynamics), and ii) unstable HBs for which d HO > 3 Å is occasionally observed (Figure 3A). Moreover, a sole H atom can form monoor multi-furcated HBs with a set of different O atoms of the zeolite framework (Figure 3B). This reflects a competition between thermally activated alkyl group motions and attractive CH…O interactions. The particular location of the H atoms with respect to the MFI silica framework would then determine the stability and strength of the HB (vide infra). Besides these geometrical and spectroscopic aspects, it is noteworthy that the covalent CH bond is polarized as deduced from the atomic partial charges (with Mulliken scheme, see Supporting Information, Tables S8), an additional criterion for the formation of HBs. [22] H atoms are positively charged on the contrary to N and C, and the differences in partial charges between H and negatively charged O are close to +0.8 confirming that the electrostatic interaction drive the CH...O bond formation. [2b,20] Therefore the OSDA-zeolite electrostatic interactions depend predominantly on the attractive interaction between H and O atoms at the organic-inorganic interface.

In the typical electron density maps of the time-averaged cluster model (Figure 4), a residual density of ca. 0.02 e -is found between H and O atoms forming HBs. This electron density explains the decrease in the global charge of TPA in interaction with the silica framework (+0.86) below its nominal charge +1. This feature, noticed previously for octadecasil and ZSM-5 zeolites containing OSDA, [9,23] suggests therefore a local stabilization of the system by charge transfers between the organic and the inorganic moieties through H-bonding. Furthermore, the three shortest and strongest CH...O bonds are spatially close to F and N, and to the plane containing these atoms (Figure 2B).

Seen along the normal axis of this plane, corresponding to the (0k0) planes in the crystal, the projections of the three HBs are located inside the triangle formed by the F atoms, each of them being between one N and one F atom (Figure 2C). This particular location appears consistent with a maximization of the electrostatic interaction between the positive charge brought by TPA and the negative charges brought by fluoride anions during the synthesis. From the potential energy density calculated at the characteristic HB critical points (3, -1), [22] the energy of the HB per O (E HB ) is estimated [24] to -1.85 kcal.mole -1 for the monofurcated bond in Figure 4A. This HB is the strongest, the shortest (2.43 Å, Figure 2C) and the most stable found by dynamics simulation (Figure 3A), and therefore one of the most able to preserve the zeolite topology during crystal growth. Its specific position, connecting a CH 2 group of TPA to the 4 MR of the [4 1 5 2 6 2 ] cage containing F, can be understood in that sense. The relative position of this CH 2 and the [4 1 5 2 6 2 ] cage is also found in optimized and non-optimized crystal structures with 2.9 < d CO < 3.7 Å (Supporting Information, Table S7b). For the other HBs E HB values are between -1.44 and -0.62 kcal.mole -1 , a typical range for weak CH…O HBs. [2b] In line with earlier H-bonding studies, [24] E HB follows an exponential law against d HO for the 28 shortest d HO that can be considered as real HBs (Supporting Information, Figure S7b). The sum of these 28 E HB values is -27 kcal.mol -1 . This energy is comparable to that calculated for the dispersion interaction of OSDA with zeolite framework (-39 kcal.mol -1 ) and represents 30% of the Coulomb interaction energy (-92 kcal.mol -1 ). In this study we have thoroughly investigated CH...O HBs in TPA-F-silicalite-1 zeolites; however short C-O distances (< 4 Å) also exist in other, very different, as-synthesized zeolites. Strikingly the normalized d CO frequency counts distributions (Figure 5) appeared rather similar for a series of assynthesized silica zeolites of contrasted topological complexity, [10] and obtained in fluoride media using organoammonium cations with different C/N ratio. This similarity can be rationalized by the favored formation of CH…O HBs, indicating that HB interactions have a widespread role in zeolite formation. AST octadecasil, [25] N,N′-diisopropylimidazolium (C/N = 4.5) for MTT, [26] hexamethonium cations (diquat, C/N = 6) for ITH-ITQ13, [27] TPA (C/N = 12) for MFI silicalite-1, [13] N,N,N-trimethyl-1-adamantammonium (C/N = 13) for STT-SSZ-23, [28] and N-benzylquinuclidinium (C/N = 14) for IFR-ITQ4. [29] OSDA-zeolite host-guest interactions are known to be governed at least by van der Waals and electrostatic interactions leading to close contacts between organic groups and the oxide network. [30] In spite of that, the existence and directional role of CH…O HBs have been largely neglected until now. Thanks to experimental data and DFT/BOMD calculations at different temperatures, we demonstrate here the ubiquitous presence of stable and dynamically unstable weak CH…O bonds in as-synthesized complex MFI zeolites, and highlighted their stabilizing role in OSDA-zeolite interactions. The identification and fine characterization of these interactions in other OSDA-zeolite framework associations would bring essential information on directing effects in zeolites. Understanding the strength and directionality of the CH…O HBs can therefore be of paramount importance for the design of the final zeolite structure and properties, such as acid site distribution and chirality.

Table of Contents

Weak CH…O hydrogen bonds between organic structure directing agents and silica framework are important host-guest interactions as demonstrated for silicalite-1 (complex topology MFI). These directional intermolecular interactions need also to be considered for various other zeolites in addition to van der Waals and Coulomb interactions.

Figure 1 .

 1 Figure 1. IR CH stretching bands of TPA-F-silicalite-1 recorded at T = 275 and 180 K.

related 19 H

 19 atoms belong to methylene or methyl groups and can form 8 monofurcated HBs (twocenter bonds with one H donor and one O acceptor) and 11 bifurcated HBs (three-center bonds with one H and two O acceptors). The histograms in Figure 2D presenting the number of d HO values below

Figure 2 .

 2 Figure 2. Time-averaged cluster: (A) general view, (B) view along zigzag channel, (C) view along the straight channel. The CH…O shortest distances are highlighted (red dashed line). The related values correspond to d HO ,

Figure 3 .

 3 Figure 3. Variations of d HO distances along the BOMD simulations for: (A) two different H and O pairs, (B) one H with three different O. For this latter case, the schemes under the curves represent the type (mono-, bi-, trifurcated or no H-bond) of H-bond interactions occurring at each time interval.

Figure 4 .

 4 Figure 4. Electron density maps of the time-averaged cluster: (A) monofurcated HB (plotted in CHO plane); (B) bifurcated HB (plotted in HOO plane). Red and blue labels correspond to atoms in plane and slightly out of plane respectively. White dots represent the critical (3,-1) points between H and O.

Figure 5 .

 5 Figure 5. Normalized frequency counts histograms of d CO distances below 5 Å retrieved from crystal structures of a series of as-synthesized zeolites. All zeolites have a silica framework and are obtained through fluoride routes. Their differences lye in the topology (see legend) and the OSDA: tetramethylammonium (C/N = 4) for

Acknowledgements

The authors thank the IRAMAN platform of the University of Montpellier for providing access to IR facility and computational resources, grant ID: A0050807369 by GENCI.