Rodrigo Méndez Rojano

Simon Mendez

Didier Lucor

Alexandre Ranc

Muriel Giansily-Blaizot

Jean-François Schved

Franck Nicoud
email: franck.nicoud@umontpellier.fr

Rodrigo Méndez

Frank Nicoud

Franck Nicoud

Kinetics of the coagulation cascade including the contact activation system: sensitivity analysis and model reduction

Keywords: Coagulation Cascade Modeling, Sensitivity Analysis, Bayesian Inference, Model reduction

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

Introduction

The use of medical devices in contact with blood has been increasing over the last years [START_REF] Sun | Antithrombotic management of patients with prosthetic heart valves: current evidence and future[END_REF][START_REF] Alderazi | Flow diverters for intracranial aneurysms[END_REF][START_REF] Indolfi | Bioresorbable vascular scaffolds -basic concepts and clinical outcome[END_REF][START_REF] Bhogal | Use of flow diverters in the treatment of unruptured saccular aneurysms of the anterior cerebral artery[END_REF][START_REF] Kirklin | Eighth annual intermacs report: Special focus on framing the impact of adverse events[END_REF]. Such devices are used to treat cardiovascular or neurovascular disorders such as coronary artery diseases, heart valve diseases or aortic/brain aneurysms. One of the main problems of bloodcontacting devices is thrombus formation which can lead to device malfunction or thromboembolism [START_REF] Mehra | The vexing problem of thrombosis in long-term mechanical circulatory support[END_REF][START_REF] Wilson | Advances in coronary stent technology: current expectations and new developments[END_REF]. Blood clotting in devices is regulated by a series of intertwined biological processes, such as protein adsorption, platelet activity, complement system and coagulation reactions, as reviewed by [START_REF] Gorbet | Biomaterial-associated thrombosis: roles of coagulation factors, complement, platelets and leukocytes[END_REF]. These mechanisms appear due to the presence of the device whose artificial material lacks the endothelial properties of the vessel. At the material surface, contact activation of factor XII (zymogen of the coagulation cascade) takes place, initiating a cascade of enzymatic reactions [START_REF] Yan | 1 -Contact activation by the intrinsic pathway of blood plasma coagulation[END_REF]) that produce thrombin, a key coagulation enzyme that activates platelets and forms fibrin. At the same time, platelet adhesion and activation can also occur at the device surface, as explained by [START_REF] Jaffer | Medical device-induced thrombosis: what causes it and how can we prevent it?[END_REF]. In the last stage of the thrombus formation process, aggregated platelets and polymerized fibrin form a stable clot that may hinder the device performance.

Computational fluid dynamics (CFD) has been used to study the flow on medical devices and evaluate the risk of thrombosis [START_REF] Yoganathan | Fluid mechanics of heart valves[END_REF][START_REF] Dumont | Comparison of the hemodynamic and thrombogenic performance of two bileaflet mechanical heart valves using a CFD/FSI model[END_REF]. A straightforward evaluation of thrombosis risk can be performed using specific flow patterns, such as high shear stress yielding to platelet activation [START_REF] Alemu | Flow-induced platelet activation and damage accumulation in a mechanical heart valve: Numerical studies[END_REF][START_REF] Shadden | Potential fluid mechanic pathways of platelet activation[END_REF] or stasis which promotes coagulation [START_REF] Biasi | Science for surgeons: Understanding pump thrombogenesis in continuous-flow left ventricular assist devices[END_REF]. Approaches focusing only on flow properties do not allow to study the biochemical processes, which play a major role in thrombus formation [START_REF] Fogelson | Fluid mechanics of blood clot formation[END_REF]. More exhaustive methods accounting for platelet activity, the coagulation cascade, clot growth and its interaction with the flow have been developed for thrombus triggered by vessel injury [START_REF] Leiderman | A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow[END_REF][START_REF] Yazdani | A general shear-dependent model for thrombus formation[END_REF]). However, only small computational domains x ∼ 100 µm can be considered due to the complexity of these models. Since a straightforward application of complicated models is not possible in the complex flow configurations linked to biomedical devices, reduced models accounting for platelet adhesion, activation (by chemical and mechanical pathways) and clot-flow interaction have been developed [START_REF] Taylor | Development of a computational model for macroscopic predictions of device-induced thrombosis[END_REF][START_REF] Wu | High fidelity computational simulation of thrombus formation in Thoratec Heart Mate II continuous flow ventricular assist device[END_REF]. Nevertheless, thrombin produced by the coagulation cascade triggered by the contact system has not been considered in a model of thrombus formation in devices.

The contact activation system has already been considered in kinetic models of the coagulation cascade [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF][START_REF] Zhu | Mathematical modeling of blood coagulation cascade: kinetics of intrinsic and extrinsic pathways in normal and deficient conditions[END_REF]. The kinetic descriptions have been used in CFD models to study thrombin formation triggered by the device wall in a typical flow configuration by [START_REF] Rojano | Introducing the pro-coagulant contact system in the numerical assessment of device-related thrombosis[END_REF]. However, a straightforward coupling with platelet-based models is challenging due to the large number of reactions present in the kinetic models and the large physical time that must be computed in a thrombus growth process. Another issue is the large uncertainty in reaction rates due to non-standardized parameter characterization or variations among different subjects [START_REF] Danforth | The impact of uncertainty in a blood coagulation model[END_REF]). In addition, according to [START_REF] Leiderman | A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow[END_REF] the use of Michaelis-Menten kinetics is questionable when enzymes and zymogens participate to several reactions. In such cases, the use of mass action descriptions of complex formation and dissociation may be more adapted to avoid any assumption on the reactions. Overall, the uncertainty related to the kinetic models may lead to poor results in thrombin production as explained by [START_REF] Belyaev | Modeling thrombosis in silico: Frontiers, challenges, unresolved problems and milestones[END_REF].

Parameter uncertainty and kinetic model reduction have been addressed for models of thrombus generation initiated by Tissue Factor exposure which is related to endothelial damage. [START_REF] Danforth | The impact of uncertainty in a blood coagulation model[END_REF] performed a sensitivity analysis on the kinetic model of [START_REF] Hockin | A model for the stoichiometric regulation of blood coagulation[END_REF] to evaluate the sensitivity of thrombin generation to the 44 reaction rates that are present in the model. The model sensitivity was assessed using the one at time (OAT) methodology with linear spaced variations between 10 to 1000% of the usual values of the reaction rates. The model outputs were evaluated at eight different instants which characterize the different phases of thrombin generation (initiation, amplification and propagation). Danforth and coworkers found that the model was especially sensitive to uncertainty in five parameters involved in the formation of the T F = V IIa complex which takes place at the initial part of the coagulation cascade. The authors suggested that improving the accuracy of the reaction rate measurement can improve the predictive capabilities of the whole kinetic model. In a posterior work, [START_REF] Danforth | Defining the boundaries of normal thrombin generation: Investigations into hemostasis[END_REF] conducted a sensitivity analysis of the initial concentrations of the eight coagulation factors involved in the model of [START_REF] Hockin | A model for the stoichiometric regulation of blood coagulation[END_REF]. They found that variations in prothrombin and T F P I accounted for 16 % and 32 % of the modifications of the model output respectively. Pairwise changes were also investigated by the authors, leading to stronger variations. The authors observed that the pair AT III with T F P I had the largest impact on the results. [START_REF] Anand | A model for the formation, growth, and lysis of clots in quiescent plasma. a comparison between the effects of antithrombin iii deficiency and protein c deficiency[END_REF] performed a similar sensitivity analysis as [START_REF] Danforth | The impact of uncertainty in a blood coagulation model[END_REF] using the model of [START_REF] Naidu | Importance of viiia inactivation in a mathematical model for the formation, growth, and lysis of clots[END_REF]. Their results showed that the generation of fibrin was most sensitive to the rate governing the production of thrombin. [START_REF] Leiderman | A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow[END_REF] recently proposed a methodology consisting of a tailored sensitivity analysis approach using the screening technique of [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] combined with a more intensive global analysis technique. The sensitivity analysis was performed for the detailed model of [START_REF] Fogelson | Blood clot formation under flow: The importance of factor xi depends strongly on platelet count[END_REF]. Three quantities of interest (QoI) related to thrombin were studied: the lag time, the maximum relative rate of generation, and the final concentration. The analysis of Link et al. (2018) considers four cases of parameter variations: variable initial concentrations levels of biochemical species, kinetic rate constants, biophysical and platelet attributes and the combined effects of varying the parameters all together. Their results show that the largest variations are observed when the combined effects of parameters variation yield to the largest effect on coagulation. Among the three studies considering only individual variations, modifications on the kinetic parameters had the largest effect on the production of thrombin.

In terms of model reduction, [START_REF] Wagenvoord | The limits of simulation of the clotting system[END_REF] showed that a model with a small number of reactions can reproduce the thrombin generation curve if it includes the basic mechanism of thrombin formation: initiation, amplification and propagation. [START_REF] Papadopoulos | A simplified mathematical model for thrombin generation[END_REF] developed a minimal model for thrombin formation based on experimental data including these essential mechanisms. Reduced-order models allow computations in complex flow configurations. For instance, [START_REF] Ngoepe | Computational modelling of clot development in patient-specific cerebral aneurysm cases[END_REF] used the reduced model of [START_REF] Wagenvoord | The limits of simulation of the clotting system[END_REF] to simulate thrombosis growth inside a cerebral aneurysm.

In the present work, a sensitivity analysis and a Bayesian inference method are used to reduce devicerelated kinetic models. First, the kinetic model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] initiated by the contact activation system is studied using a global screening technique that allows to identify the most sensitive parameters. Then, a reduced chemical model triggered by contact activation of factor XII is proposed. The optimal parameters of the reduced model are obtained using a data-driven Bayesian statistical framework to find the optimal parameters used in the reduced model.

Methods

Plasma Samples and TGA

Thrombin generation assay (TGA) was used to measure the evolution of thrombin concentration in time after coagulation was triggered by the contact activation system. Developed by [START_REF] Hemker | Calibrated automated thrombin generation measurement in clotting plasma[END_REF], the Calibrated Automated Thrombogram (CAT) is a test mainly used in hemostasis research to study the hemostatic profile. This method is time-consuming and suitable only for a small number of samples when compared to routine clinical tests which assess the first traces of thrombin (after coagulation is triggered) about 12 seconds for Quick Time (QT) or 32 seconds for Activated Partial Thromboplastin Time (APTT) for a normal plasma. In contrast, TGA assesses thrombin generation until 60 minutes and allows a better characterization of the coagulation cascade throughout all its phases: initiation amplification and propagation. Any traces of thrombin will cleave to a specific substrate Z-Gly-Gly-Arg-AMC into a flurogenic form, detected by the thrombogram.

Samples: Each biological assays was performed with the Pooled Normal Plasma (PNP R ) [Cryopep, Montpellier France] which consists of and industrial pooled citrated platelet poor plasmas (PPP) obtained from at least 20 healthy donors to avoid inter-individual physiological variations of clotting factor levels. Each donor has been pretested for QT, APTT, factor II, V, VII, VIII, IX, X, XI, XII levels, which all were between normal range. The PNP was diluted with prothrombin immuno-deficient plasma whose qualified activity was lower than 1% [Siemens Healthcare, Erlangen, Germany] in order to obtain different final prothrombin concentrations. Silica was used to trigger the coagulation cascade, mixed with rabbit cephalin in STA-PPT A [Stago, Asnieres-sur-Seine, France].

Calibrated automated measurement of thrombin generation: Thrombin generation was determined in PPP using Fluoroscan Ascent (Flucakit, Thrombinoscope, Synapse BV, Maastrich, The Netherlands) according to the method described by [START_REF] Hemker | Calibrated automated thrombin generation measurement in clotting plasma[END_REF]. 80 µL of plasma were mixed with 20 µL of STA-PPT A and 20 µL of fluorescent reagent, FlucaKit. This reagent contains calcium chloride (necessary to trigger the coagulation cascade) and a mixture of fluorogenic substrate (Fluo-Substrate) and FluoBuffer. Thus, fluorescence intensity was detected at wavelengths of 390 nm (excitation filter) and 460 nm (emission filter), every 20 seconds. Each individual sample is analyzed with a thrombin calibrator as reference for a stable thrombin activity of approximately 600 nM. The calibrator enables the conversion of the fluorescence signal into thrombin concentration. The signal is treated to correct inner filtering effects, substrate consumption and abnormal plasma color. Analyses were conducted, on Immulon 2HB round-bottom 96-well plates (Stago -Asnières-sur-Seine, France).

Detailed coagulation model

Sensitivity analysis is performed on the model introduced by [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]. The model includes 37 reactions and 63 reaction rates and features:

-Initiation by the contact activation system considering factor XII auto-activation, reciprocal activation (by Kallikrein) and auto-hydrolysis (by factor XII a ), -Extrinsic and common pathways following the model of [START_REF] Hockin | A model for the stoichiometric regulation of blood coagulation[END_REF], -Inhibition including Antithrombin (AT III), Kallikrein inhibitor K inh and C 1 inhibitor C1 inh .

The TGA experimental data generated in our laboratory was use as a test case for the model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] with some modifications related to

# Reaction Forward M -1 s -1 Reverse s -1 Forward s -1 1 Xa + V II → Xa + V IIa k 1 = 1.3 × 10 7 2 IIa + V II → IIa + V IIa k 2 = 2.3 × 10 4 3 II + Xa → IIa + Xa k 3 = 7.5 × 10 3 4 IIa + V III → IIa + V IIIa k 4 = 2.0 × 10 7 5 V IIIa + IXa ↔ IXa = V IIIa k 5 = 1.0 × 10 7 k 6 = 5.0 × 10 -3 6 IXa = V IIIa + X ↔ IXa = V IIIa = X → IXa = V IIIa + Xa k 7 = 1.0 × 10 8 k 8 = 1.0 × 10 -3 k 9 = 8.2 7 V IIIa ↔ V IIIa 1 • L + V IIIa 2 k 10 = 6.0 × 10 -3 k 11 = 2.2 × 10 4 8 IXa = V IIIa = X → V IIIa 1 • L + V IIIa 2 + X + IXa k 12 = 1.0 × 10 -3 9 IXa = V IIIa → V IIIa 1 • L + V IIIa 2 + IXa k 13 = 1.0 × 10 -3 10 IIa + V → IIa + V a k 14 = 2.0 × 10 7 11 Xa + V a ↔ Xa = V a k 15 = 4.0 × 10 8 k 16 = 0.2 12 Xa = V a + II ↔ Xa = V a = II → Xa = V a + mIIa k 17 = 1.0 × 10 8 k 18 = 103 k 19 = 63.5 13 Xa = V a + mIIa → Xa = V a + IIa k 20 = 1.5 × 10 7 14 Xa + T F P I → Xa = T F P I k 21 = 9.0 × 10 5 k 22 = 3.6 × 10 -4 15 Xa + AT III → Xa = AT III k 23 = 1.5 × 10 3 16 mIIa + AT III → mIIa = AT III k 24 = 7.1 × 10 3 17 IXa + AT III → IXa = AT III k 25 = 4.9 × 10 2 18 IIa + AT III → IIa = AT III k 26 = 7.1 × 10 3 19 BocV P RM CA + IIa → BocV P RM CA = IIa k 27 = 1.0 × 10 8 k 28 = 6.1 × 10 3 k 29 = 53.8 20 XII → XIIa k 30 = 5.0 × 10 -3 21 XIIa + XII ↔ XIIa = XII → XIIa + XIIa k 31 = 1 × 10 8 k 32 = 750 k 33 = 3.3 × 10 -2 22 XIIa + P K ↔ XIIa = P K → XIIa + K k 34 = 1 × 10 8 k 35 = 3.6 × 10 3 k 36 = 40 23 XII + K ↔ XII = K → XIIa + K k 37 = 1 × 10 8 k 38 = 45.3 k 39 = 5.7 24 P K + K → K + K k 40 = 2.7 × 10 4 25 K → Kinh k 41 = 1.1 × 10 -2 26 XIIa + C1inh → XIIa = C1inh k 42 = 3.6 × 10 3 27 XIIa + AT III → XIIa = AT III k 43 = 21.6 28 XI + IIa ↔ XI = IIa → XIa + IIa k 44 = 1 × 10 8 k 45 = 5 k 46 = 1.3 × 10 -4 29 XIIa + XI ↔ XIIa = XI → XIIa + XIa k 47 = 1.0 × 10 9 k 48 = 200 k 49 = 5.7 × 10 -3 30 XIa + XI → XIa + XIa k 50 = 3.19 × 10 6 31 XIa + AT III → XIa = AT III k 51 = 3.2 × 10 2 32 XIa + C1inh → XIa = C1inh k 52 = 1.8 × 10 3 33 XIa + α 1 -AT → XIa = α 1 -AT k 53 = 1.0 × 10 2 34 XIa + α 2 -AP I → XIa = α 2 -AP k 54 = 4.3 × 10 3 35 XIa + IX ↔ XIa = IX → XIa + IXa k 55 = 1.0 × 10 8 k 56 = 41.0 k 57 = 7.7 36 IXa + X ↔ IXa = X → IXa + Xa k 58 = 1.0 × 10 8 k 59 = 0.64 k 60 = 7.0 × 10 -4 37 Xa + V III ↔ Xa = V III → Xa + V IIIa k 61 = 1.0 × 10 8 k 62 = 2.1 k 63 = 0.023
Table 1 Kinetic parameters and coagulation reactions used in the simulations. The model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] was adapted to the experimental conditions detailed in Section 2.1.

our experimental setup. The inhibition by corn trypsin inhibitor was suppressed since it was not used. The initial trigger reaction activation of factor XII (reaction 20) was modified, the original value of Chatterjee et al. ( 2010) being fitted to their specific experimental framework. Fibrin related reactions were not included in the model to reduce the computational cost of the model when performing the sensitivity analysis. The reactions related to factor X and factor IX activation by factor V II a were deleted because those reactions referred to the non-physiological recombinant factor V II a [START_REF] Komiyama | Proteolytic activation of human factors ix and x by recombinant human factor viia: effects of calcium, phospholipids, and tissue factor[END_REF]) (these reactions are numbered 48 and 49 in [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]). The prediction of the resulting model is plotted in Fig. 1 in dash-dotted line. It was observed that the model systematically overestimates the lag time. In order to improve the baseline comparison between experiments and model predictions, the kinetic rates k 47 and k 49 were augmented by a factor of 10 with respect to the values reported by [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]. The result is plotted in Fig. 1 as solid line and compares well with the experimental data. All the reactions and the baseline values of the reaction rates are reported in Table 1. Note that these modifications drastically improve the numerical-experimental comparison in the thrombin generation test cases considered in this study.

A system of ordinary differential equations (ODE) can be written from the biochemical reactions of the kinetic model (Table 1). The ODE system writes:

dC i dt = R i (1)
where C i is the concentration of species i and R i is the source term which is obtained by applying the law of mass action to the coagulation reactions of Table 1. As an example the source term of Kininogen K is devel-oped as:

dC K dt = -k 41 C K -k 37 C XII C K + (k 38 + k 39 )C XII=K + 2k 40 C P K C K + k 36 C XIIa=P K (2)
The solution of the system is computed using an integration Runge-Kutta scheme of the 4th order. This numerical scheme was implemented in the in-house CFD solver YALES2BIO 1 . YALES2BIO has already been used in coagulation problems by [START_REF] Rojano | Introducing the pro-coagulant contact system in the numerical assessment of device-related thrombosis[END_REF] and other cardiovascular flows in micro [START_REF] Sigüenza | How should the optical tweezers experiment be used to characterize the red blood cell membrane mechanics?[END_REF][START_REF] Lanotte | Red cells' dynamic morphologies govern blood shear thinning under microcirculatory flow conditions[END_REF]) and macro [START_REF] Chnafa | Image-based simulations show important flow fluctuations in a normal left ventricle: What could be the implications?[END_REF][START_REF] Zmijanovic | About the numerical robustness of biomedical benchmark cases: Interlaboratory FDA's idealized medical device[END_REF][START_REF] Sigüenza | Fluid-structure interaction of a pulsatile flow with an aortic valve model: A combined experimental and numerical study[END_REF]) scale applications.

Table 2 shows the initial concentrations used in the numerical simulations which correspond to the physiological concentrations of coagulation factors after sample preparation.

In addition to the concentrations of species, which are patient-specific, the coagulation model is governed by 63 reaction rates (see Table 1). The reaction rates are derived from reactions observed under experimental conditions of saturating phospholipid concentration. Most of the kinetics models are build upon reaction rate values reported in literature [START_REF] Jones | A model for the tissue factor pathway to thrombin. ii. a mathematical simulation[END_REF][START_REF] Zhu | Mathematical modeling of blood coagulation cascade: kinetics of intrinsic and extrinsic pathways in normal and deficient conditions[END_REF]; [START_REF] Zarnitsina | A mathematical model for the spatio-temporal dynamics of intrinsic pathway of blood coagulation. i. the model description[END_REF]. Additional fitting is in general required due to non-physiological conditions in which some parameters are derived or extensively modified proteins which are used to obtain the parameter value [START_REF] Hockin | A model for the stoichiometric regulation of blood coagulation[END_REF]. For this reason, reported values in different kinetic models sometimes present large variation. For instance, the reaction rate of factor XI activation can be found in literature with a difference up to 5 orders of magnitude as highlighted by [START_REF] Belyaev | Modeling thrombosis in silico: Frontiers, challenges, unresolved problems and milestones[END_REF].

Figure 1 shows the numerical and experimental data. The numerical set up here presented was used as a nominal reference for the sensitivity analysis. A reasonable agreement between the experimental TGA and the kinetic model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] was obtained.

Global Sensitivity Analysis by Morris screening

The objective of the present sensitivity analysis is not to precisely quantify the influence of the kinetic parameters on coagulation as in [START_REF] Danforth | The impact of uncertainty in a blood coagulation model[END_REF] kinetic rates. This is a necessary step in order to gain more insight into the model and use this knowledge to build up the reduced kinetic model of Section 3.2. The screening method of [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] allows a fast exploration of a model through the discretization of the input parameters. This type of method is well suited to models with many parameters and a good compromise between accuracy and efficacy has been reported [START_REF] Iooss | Uncertainty Management in Simulation-Optimization of Complex Systems: Algorithms and Applications[END_REF]. The Morris method was also recently used by [START_REF] Leiderman | A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow[END_REF], who highlighted its ability to rank the most sensitive parameters of a complex thrombosis model. The general idea of the Morris analysis is to calculate a number of elementary effects for each input parameters and to compute basic statistics to identify the most sensitive parameters. In order to compute the elementary effects a discrete grid is built on the input variables and explored with an efficient sampling technique. In the current work, the python module SALIB 2 that follows the optimized version of the Morris method by [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF] was used. We considered only the case of variations on the kinetic rates as they produced the largest change in the production of thrombin as observed in [START_REF] Leiderman | A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow[END_REF]. The practical steps that were performed to conduct the analysis using the model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] are the following:

1. Choice of discrete grid and sampling: a discrete grid on the p = 63 input variables is built with l probability levels uniformly distributed between [0, 1]. In this analysis 10 levels of discretization were used.

Inside the input domain a number of trajectories r has to be determined in order to cover the input parameters space in an efficient manner. In the original method of [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF], the trajectories starting point is chosen randomly, then consecutive one-at-a-time increments are performed in the discrete grid with a random direction. To perform an optimal coverage of the input space, the optimized strategy of [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF] was preferred in this analysis. In the optimal strategy, the dispersion of the starting points is maximized through the generation of several Morris trajectories (in this case r max = 500). Then the most spread optimal trajectories are kept, for this case r = 10. The total number of experiments is then defined by M = r(p + 1) thus, for the current model, M = 640 TGA simulations were required. 2. Mapping of parameters to actual distribution: once the sampling is performed a Design of Experiment (DoE) matrix with parameters between 0 and 1 is obtained. A quantile function is applied to transform the parameters to their actual values. This is a common practice in several screening strategies [START_REF] Saltelli | Sensitivity analysis practice: A guide to scientific models[END_REF]. In our analysis, a loguniform probability distribution is used for all the input parameters. This choice was motivated by the large parameter variations reported in the literature [START_REF] Hemker | Is there value in kinetic modeling of thrombin generation? no (unless[END_REF]). In addition, the actual distribution for each reaction rate is not available in literature. The quantile function used writes

Q(z i ) = c inf,i c sup,i c inf,i zi (3)
where [c inf,i ; c sup,i ] is the distribution support and z i is the sampled value of the i parameter (between [0,1]). To guarantee that the mean µ i of the sampled parameters coincides with the nominal input value k i , the superior and the inferior support values must be estimated. To do so, the support ratio β = csup,i c inf,i 2 http://salib.readthedocs.io/en/latest/api.html is introduced in the log-uniform mean expression

µ i = (c sup,i -c inf,i )/(log(c sup,i ) -log(c inf,i )). Im- posing µ i = k i then leads to c inf,i = kilog(β)
β-1 and c sup,i = βc inf,i . In the current analysis, β was set equal to 100 which translates in parameters variation of about 4.5% -450% of each parameter nominal values. Two trajectories in the (k 3 , k 7 ) plane are displayed in Figure 2 to illustrate the method.

E i j = f (X i j + δ i ) -f (X i j )) δ i (4)
where X is the input parameter, j is the index of the input parameter, i stands for the repetition, δ i is the variation value and the function f represents any of the five QoI from the thrombin generation curve. Two sensitivity coefficients are calculated [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF]):

-µ * j = 1 r r i=1 | E i j |
is the mean of the absolute value of the elementary effects

-σ j = 1 r r i=1 E i j -1 r r i=1 E i j 2
is the standard deviation of the elementary effects.

From the sensitivity indices distribution a quick identification of three different groups of parameters can be done depending on their impact on the QoI:

parameters with little or no effect (low µ * j and σ j ), parameters with a strong linear dependency or an additive effect (large µ * j ), parameters with non-linear and/or interaction effects (large σ j ), With this classification the most sensitive parameters can be identified. In some cases, the sensitivity analysis can identify predominant reactions and contribute to reduce the complexity of the models, as discussed in Section 4.

2.4 A reduced model of coagulation initiated by contact activation [START_REF] Wagenvoord | The limits of simulation of the clotting system[END_REF] proposed a minimal reaction mechanism of thrombin formation triggered by factor V II a . The reduced model included reactions of the extrinsic and common pathways and can thus be applied in cases were T F is the initial mechanism. In this model, 14 parameters including 9 kinetic constants and 5 initial concentrations of coagulation factors must be optimized in order to correctly reproduce experimental trends of thrombin formation. It is worth mentioning that in reduced kinetic models, the species maybe fictive species and the optimal initial concentrations of factors are not representative of the actual concentrations in blood plasma. However, these models have been successful to reproduce thrombin formation triggered by the extrinsic pathway. In the current work a reduced model for Fig. 3 Variables used in the sensibility analysis:

t lag , max(II a ), m 1 , m 2 , ETP. XII XII a II II a V V a AT III AT III = II a k 1 K m2 k cat2 K m3 k cat3 K m4 k cat4 k 5
Fig. 4 Reduced model for thrombin generation triggered by the contact activation system. The notation K mi , k cati is used for reaction following MichaelisMenten kinetics while k i are first order reactions.

thrombin production initiated by contact activation is proposed. The minimal reaction mechanisms consist of 5 reactions that are triggered by activation of factor XII, the reactions are listed on Table . 3. The model includes 8 reaction rates and 8 chemical species. The amplification and propagation phases are performed by the auto-activation loop through the activation of factor V a , this mechanism also appearing in the model of [START_REF] Wagenvoord | The limits of simulation of the clotting system[END_REF]. Finally, inhibition of thrombin is due to the activity of AT III. Figure 4 shows a schematic of the model with the corresponding kinetic parameters. The kinetic model Ċ = M(C, t) can be

# Reaction K m M k cat s -1 k s -1 1 XII → XII a k 1 2 XII a + II → II a + XII a K m2 k cat2 3 II a + V → II a + V a K m3 k cat3 4 V a + II → V a + II a K m4 k cat4 5 II a + AT III → II a = AT III k 5
Table 3 Reduced kinetic model for the coagulation cascade showed in Fig. 4. Reactions 1 and 5 are first order reaction, reactions 2, 3 and 4 follow MichaelisMenten kinetics.

represented by a system of ordinary equations that can be solved using a classical numerical scheme. The ODE system writes:

                             ĊXII ĊXIIa ĊII ĊIIa ĊV ĊV a ĊAT III ĊAT III=IIa                              =                               -k 1 C XII k 1 C XII k cat4 C Va C II K m4 + C II - k cat2 C XIIa C II K m2 + C II -k 5 C IIa C AT III k cat2 C XIIa C II K m2 + C II + k cat4 C Va C II K m4 + C II -k 5 C IIa C AT III - k cat3 C IIa C V K m3 + C V k cat3 C IIa C V K m3 + C V -k 5 C IIa C AT III k 5 C IIa C AT III                              
Parameter values should be optimized to compute thrombin formation with a reasonable agreement with reference data. Here we use the experimental data generated by TGA triggered by contact activation.

Bayesian parameter inference

The proposed reduced-order coagulation model depends on a set of parameters which need to be calibrated. The idea behind this parametric calibration is to adjust the unknown parameters in order to lower the discrepancy between the numerical prediction and some indirect and imperfect (experimental) observations of the system. This is a computational burden because the model is nonlinear and brings in unknown relatively high-dimensional parameter spaces. Model inversion in the presence of measurement errors must typically take advantage of some type of regularization (e.g., Tikhonov regularization) in order to recover the existence and uniqueness of solutions or a robust optimization method [START_REF] Ashyraliyev | Systems biology: parameter estimation for biochemical models[END_REF]. A potentially more natural setting is the Bayesian statistics. The Bayesian statistical framework can be explained as a systematic use of probability dedicated to decisionmaking when using a model with uncertain parameters.

It is advantageous because it allows the specification of a prior distribution which expresses probabilistically what is known about the parameters before observing the data. In this study, a Bayesian inference technique is used to obtain the optimal parameters involved in the model that allow to capture the thrombin formation curve. The parameters to be inferred are the 8 reactions rate and the initial concentrations of species XII, V and AT III, i.e. a total of 11 quantities. This approach is used in inverse problems and relies on posterior sampling techniques. Thus, a prior distribution of the parameters to be calibrated is necessary, so that the estimation process delivers a probabilistic characterization of the parameters. In this work, a Markov-Chain Monte-Carlo algorithm is used considering the input parameters as random quantities that are sampled according to the posterior distribution.

A significant amount of literature exists on the Bayesian inference method. In the following, the work of Xiu (2010); [START_REF] Birolleau | Adaptive bayesian inference for discontinuous inverse problems, application to hyperbolic conservation laws[END_REF]; [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF] is referenced to recall the main steps of the parameter inference. In the Bayesian approach, a vector of unknown parameters is considered

k = (k 1 , K m2 , k cat2 , K m3 , k cat3 , K m4 , k cat4 , k 5 , XII init , V init , AT III init ).
(5)

The parameters are treated as random variables with a prior density distribution π prior (k). In the current inference procedure, independent log-uniform distributions were assumed as prior distributions for all the parameters. Observations of the QoI are necessary to perform the inference process. In this application, the observed variables d ∈ R n d are the experimental measurements of the thrombin generation assay or the numerical results of the detailed kinetic model at n t discrete time instants C IIa (n t ).

The Bayesian formula, involving conditional probabilities, can be applied as in [START_REF] Birolleau | Adaptive bayesian inference for discontinuous inverse problems, application to hyperbolic conservation laws[END_REF] to obtain the parameters posterior distribution considering the reduced model, the prior parameters distribu-tion and the experimental observations

π post (k|d) ∝ π (d|k)π prior (k), (6)
where π (d|k) is the likelihood function that combines the experimental data with the forward model and π post (k|d) the desired posterior density of the parameters. The initial concentration of prothrombin is not inferred, thus the physiological condition of pool plasma is used in the reduced kinetic model, II = 933.0 nM. Following [START_REF] Xiu | Numerical Methods for Stochastic Computations[END_REF] assuming additive measurement noise ε and mutually independent random variables, one can write:

d = G(C) + ε = G (M(k, t)) + ε, (7)
where

M : R n k × R nt → R n C ×nt is the deterministic forward model, G : R n C ×nt → R n d is an observation
operator that relates the model solution C with the reference discrete concentration of thrombin d and π ε is the prescribed noise distribution of ε. The observation operator G can be any operation that transforms the model outputs into the QoI used in the inference process. In the current case the operator G only selects the thrombin concentration C IIa from the output concentration vector C of the forward model M. More specifically, in our case, the likelihood function can be expressed as follows:

π (d|k) = n d i=1 ε (d i -G i (M(k, t))) . ( 8
)
where the index i is used to refer to the experimental and numerical observations at the same n d time instants. The likelihood function contains a stochastic source term that must encompass the response of the deterministic kinetic model over the support of π post (k|d). Since the posterior distribution does not have in this case an analytic closed expression, posterior moments, expectations or maximum a posteriori values must be estimated via sampling methods such as Markov chain Monte Carlo (MCMC) as explained in [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF]. The posterior sampling of the parameters is handled thanks to a standard Metropolis-Hastings scheme [START_REF] Hastings | Sources of circuit thrombosis in pediatric extracorporeal membrane oxygenation[END_REF]). In our case, it requires many model realizations (about twenty thousand), which is not penalizing as the reduced-order model is not computationally demanding.

Results

Sensitivity Analysis of the full coagulation model

The first aspect of interest is the robustness of the detailed coagulation model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] ini-tiated by contact activation. Figure 5 shows the results of the sensitivity analysis using the screening Morris method for each quantity of interest. The basic statistics of the elementary effects σ j and µ * j are computed using the thrombin concentration of the 640 simulations. An arbitrary threshold of 30% of the largest value of both sensitivity indices is used to keep apart the sensitive parameters.

The important parameters identified for the time lag t lag (Fig. 5a) belong to the intrinsic pathway; activation of factors IX (reaction 35) and XI (reaction 29) have large impacts on t lag . The parameters relevant to max(II a ) are reactions rates involved in the activation of factor II (Fig. 5b). In addition, thrombin inhibition due to AT III presents a large influence on the maximum value of thrombin. The ascending slope m1 (Fig. 5c) shows a large influence to reaction rates involved in reaction 12 which is formation of mII a by the Prothrombinase complex. Furthermore, reaction rates involved in the intrinsic pathway, activation of factors, XI a , IX, X and V III, showed to be important for m 1 for instance. The descending slope of thrombin evolution m 2 (Fig. 5d) is mostly influenced by the thrombin inhibition activity of AT III and by the production of mII A and II a in reactions 12 and 13 in which kinetic rates 18 and 20 participate.

Finally, the results on ETP (Fig. 5e) show that the most influential parameters are involved in the inhibition activity of AT III by means of complex formation with II a and mII a . Reaction rate k 20 which is involved in the conversion of mII a into II a is also relevant to ETP. The histogram shows that the conversion of mII a into II a determined by k 20 is relevant to all the QoI. In a similar way, inhibition activity by AT III due to parameter k 26 is important for 4 QoI, excluding only t lag . Six reaction rates are sensitive to at least three QoI, the majority of these parameters are involved directly in the production of II a with the exclusion of parameter k 60 which is involved in the activation of factor X by IX a .

The parameters identified by the sensitivity analysis are consistent with the nature of the QoI. It is important to note that these results were verified using a very different strategy for global sensitivity analysis. In this case, influential parameters were identified thanks to their first-order Sobol' indices [START_REF] Sobol | Global sensitivity indices for rather complex mathematical models can be efficiently computed by Monte Carlo (or quasi-Monte Carlo) methods. these indices are used for estimating the influence of individual variables or groups of variables on the model output[END_REF]). Sobol' indices are a variance-based diagnostic method of global sensitivity analysis that can identify the parameters that cause large variations on the output of a given model. While those indices may be evaluated with Monte-Carlo sampling type method, they are here approximated from a surrogate polynomial approximation of each QoI. These polynomial approximation were built from level-3 Stroud quadrature sampling via numerical projections on Legendre first-order orthogonal polynomials of k [START_REF] Xiu | Numerical Methods for Stochastic Computations[END_REF]. In this case, only 126 simulations were performed. Figure 6 shows the Sobol' coefficients of the 63 parameters for t lag . The parameters marked with a cross are the parameters with large variance. In fact, a remarkable consistency with the results from the Morris method is observed for all the QoI (results for the 4 remaining QoI are not shown).

0

Reduced model for thrombin generation

Whatever the confidence we have in a coagulation model, it can only be used in CFD if its computational cost remains tractable. As a result, model reduction is particularly useful to build coagulation models with a limited number of species. Of course, such models involve new parameters whose values need to be determined using reference data sets.

Eight thrombin generation assays initiated by contact activation were conducted on the same PNP sample to obtain the data sets relevant to thrombin. The mean of these thrombin formation curves was used to perform the Bayesian inference assuming a 4 % of observation uncertainty. Optimal parameters obtained after 20, 000 MCMC (Metropolis-Hastings sampling algorithm) iterations are listed in Table 4 Figure 7 shows the experimental data and the numerical results using the Bayesian optimal parameters in the reduced coagulation model. Numerical results using the Bayesian optimal parameters aligned well with TGA experimental data. To evaluate the predictive capabilities of the proposed reduced model, the initial condition of factor II was varied. Figure 8 shows both the numerical and TGA results for three prothrombin concentrations. The numerical results reproduce fairly well the trends of thrombin formation, demonstrating the good predictive capabilities of the reduced model using the Bayesian optimal parameters. The three cases correspond to the physiological condition C II = 933.0 nM, 50% and 15% of the physiological values, respectively. It can be noticed that in both cases at low prothrombin concentrations are shifted with longer lag times than the TGA data. In addition, the amplitude of the II a curves is, in both cases, smaller than the experimental results. The inference process was also repeated with a different set of reference data using the numerical data from the detailed model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF] assuming an uncertainty on thrombin concentration of 4 %. Table 5 shows the values of the inferred parameters. This approach can be applied to more general cases in which the experimental conditions differ from the TGA samples previously presented.

Figure 9 shows the thrombin generation curves using the inferred parameters of 2010) and numerical data using the reduced model with optimal parameters of Table 5 for a range of factor II initial concentrations, 100% =physiological PPP concentration, 50%, 15%)

Discussion

In this work, two kinetic models of the coagulation cascade initiated by contact activation with a different level of complexity are used: a full reference model from the literature and a reduced model introduced here for the first time. First, the [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] method for sensitivity analysis is applied to the existing detailed model of the coagulation cascade initiated by the contact ac-tivation system. This screening technique allowed the identification of the parameters that have a significant impact on the model QoI's. The kinetic rates that had the largest impact on the time before the reactions take up t lag are involved in the intrinsic coagulation pathway (k 57 , k 55 , k 60 , k 56 , k 48 , k 49 ). Interestingly, the reaction rate k 30 which is implicated in the contact activation of factor XII does not have a significant effect on t lag , thus suggesting that its role is only to trigger the cascade of reactions. This results are consistent with the numerical results of Méndez [START_REF] Rojano | Introducing the pro-coagulant contact system in the numerical assessment of device-related thrombosis[END_REF] in which a parametric study on the activation rate of factor XII a showed a limited impact on thrombin formation. This effect may explain the difficulties of developing a full non-thrombotic material since reducing the contact activation rate does not help to reduce the thrombin concentrations levels.

The ascending slope m 1 is particularly sensitive to k 50 which is involved in the auto-activation of factor XI. In addition parameters k 61 and k 63 participating in the formation of thrombin and activation of factor V III showed large impact on this QoI. For the other quantities max(II a ), m 2 and ETP, the reaction rates that were identified correspond to reactions that are directly involved in the production of thrombin or the inhibition activity of AT III (k 26 , k 20 , k 23 , k 24 , k 18 , k 19 ).

To verify the results, another sensitivity analysis relying in the Stroud sampling quadrature and the Sobol' indices was performed and the results were qualitatively confirmed for all the QoI's. In the histogram showed in Fig. 5, parameter k 20 was identified to be important for all the QoI. Reaction rate k 26 is sensitive to four QoI except for t lag . These findings align well with the sensitivity study of [START_REF] Danforth | The impact of uncertainty in a blood coagulation model[END_REF] on the coagulation model triggered by T F = V II a . Danforth and coworkers found that the thrombin formation levels were sensitive to variations on reaction rates k 26 and k 20 . Many reaction rates identified by Danforth participate in the initiation phase of the extrinsic pathway. In the current work reactions rates k 57 , k 56 and k 55 participating in the intrinsic system have a considerable impact on the QoI's t lag and m 1 . Finally, six reactions rates were found to be important to three output variables, k 18 , k 19 , k 24 which are involved in the production and inhibition of mII a and k 60 which participates in the activation of factor X by IX a . A better experimental characterization of these six reactions may lead to an improvement in the robustness of the model. We emphasize that other uncertainty quantification techniques more accurate than Morris method also exist for a moderate number of random parameters (?). However, the parameters identified in the current work can be used to perform a more accurate sensitivity analysis on the most important parameters selected thanks to the Morris method, as done for instance by [START_REF] Leiderman | A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow[END_REF].

As pointed out by [START_REF] Ngoepe | Thrombosis in cerebral aneurysms and the computational modeling thereof: A review[END_REF] thrombosis models should be sophisticated enough to provide precious clinical information with a reasonable computational cost. To tackle this challenge in the context of deviced-associated thrombosis, a reduced model for thrombin generation triggered by contact activation of factor XII is proposed. This model allows for the characterization of the patient thrombin formation profile with a minimal reaction mechanisms. The reduced model was constructed based on previous reduction approaches [START_REF] Wagenvoord | The limits of simulation of the clotting system[END_REF] and by substituting the extrinsic pathway with the contact activation system. In fact, the sensitivity results also guided the construction of the reduced model. For instance, the sensitivity analysis showed that formation of thrombin and meizothrombin by factor V a had an important role by means of parameters k 18 , k 19 and k 20 . In addition, inhibition activity of AT III showed to be relevant through parameters k 26 and k 24 . Thus, these two mechanisms were incorporated with the minimum number of reactions in the reduced model. Finally, in order to compute thrombin evolution, the 11 parameters involved in the model (kinetic rates and initial constants) have to be calibrated. A Bayesian inference framework with a MCMC sampling method was used to infer the optimal parameters using the concentration of thrombin in time as reference data from both the TGA and the detailed model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]. The set of parameters found by the Bayesian inference method could reproduce the experimental TGA trend when the initial concentration of factor II was decreased to 50 and 15 % from its physiological value. A supplemental set of parameter values was inferred using the detailed kinetic model of Chatterjee et al. The set of parameters aligned well with the trend presented by Chatterjee's model; this set of parameters may be used in general purpose CFD computations in which the conditions of the case are not the same as in our experimental configuration. Both sets of optimized parameters showed to follow a good trend with respect to TGA and numerical data showing equivalent predictive capabilities. Using such a reduced kinetic model in hemodynamic computations lowers the computational cost of chemical species transport by almost a factor of seven.

The results contribute to improve the predicting capabilities of the model of [START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]. To the best of our knowledge, this is the first time that a minimal reaction mechanism triggered by the contact activation system is proposed. The reduced model combined with proper boundary conditions [START_REF] Rojano | Introducing the pro-coagulant contact system in the numerical assessment of device-related thrombosis[END_REF]) may be particularly useful in situations in which the coagulation reactions are triggered by contact activation of factor XII.

Fig. 2

 2 Fig.22D parameter grid, the levels are mapped using the log uniform distribution. Two trajectories are shown as an ilustration of the OAT random advancement.

Fig. 5

 5 Fig.5Each QoI statistics ( σ j vs µ * j ) are evaluated through the Morris analysis. The rectangle represents the threshold retention zone, the parameter with effects larger than 0.3 × max(σ j ) and 0.3 × max(µ * j ) are considered relevant. A histogram showing how many QoI are sensitive to each parameter is also shown (bottom right).

Fig.

  Fig.Mean of the eight TGA data sets for physiological data and numerical thrombin formation using the reduced model (Fig.4) with optimal parameters from Bayesian inference. The error bars represent the variance of the TGA data sets. The optimal parameters were obtained after 20000 MCMC iterations.

Fig.

  Fig. Numerical results of the detailed model of Chatterjee et al. (2010) and numerical data using the reduced model with optimal parameters of Table5for a range of factor II initial concentrations, 100% =physiological PPP concentration, 50%, 15%)

  [START_REF] Leiderman | A local and global sensitivity analysis of a mathematical model of coagulation and platelet deposition under flow[END_REF], but rather to identify the most influencing 1 http://www.math.univ-montp2.fr/∼yales2bio/

			Factor	Initial Condition [nM]
			V II		6.67
			V IIa		0.667
			X		106.67
			IX		6.0
			II		933.0
			V III		0.4667
			V		13.33
			T F P I		1.667
			AT III		2267.0
			XII		226.7
			P K		300.0
			C1 inh XI		1667.0 206.7
			α 1 -AT	30000.0
			α 2 -AP	667.0
			BocVPRMCA	6670.0
	Table 2 Factors concentrations based on TGA (concen-
	trations are calculated after sample mixing with fluorescent
	reagent and using the physiological values reported by Chat-
	terjee et al. (2010)) used as initial conditions in the numerical
	simulations.
			•10 -7
			6
		II a (M)	2 4
			0
			0	100	200	300	400
				time (s)
	Fig. 1 Experimental results of thrombin generation assay
	( ) and numerical simulation (	) using the kinetic scheme
	of Table 1 with the initial conditions of Table 2. The nu-
	merical result without modifying k 47 and k 49 is also shown
	(	).

  .

	#	k	optimal
	1	k 1	6.4 × 10 -3
	2 3 4 5	k cat 2 k m 2 k cat 3 k m 3	8.21 8.95 × 10 -9 3.5 × 10 -3 2.0 × 10 -9
	6 7 8	k cat 4 k m 4 k 5	4.98 8.25 × 10 -7 7.79 × 10 3
	9	C XII	9.36 × 10 -11
	11	C V	6.2 × 10 -9
	12 C AT III	1.665 × 10 -6

Table 4

 4 Optimal parameters obtained from Bayesian inference; (here taken as the statistical mode values) the initial concentration of prothrombin C II = 933.0 nM is not infered.

Table 5 .

 5 A lag time between the model of Chatterjee et al. and the reduced model is observed for the case of 15% prothrombin ini-

	#	k	optimal
	1	k 1	1.2 × 10 -3
	2 3 4 5	k cat 2 k m 2 k cat 3 k m 3	0.28 3.72 × 10 -7 7.7 × 10 -3 6.4 × 10 -10
	6 7 8	k cat 4 k m 4 k 5	7.01 9.54 × 10 -7 9.88 × 10 3
	9	C XII	7.91 × 10 -11
	11	C V	1.8 × 10 -8
	12 C AT III	1.977 × 10 -6

Table 5

 5 Optimal parameters obtained from Bayesian inference using the numerical data of the detailed model of[START_REF] Chatterjee | Systems biology of coagulation initiation: Kinetics of thrombin generation in resting and activated human blood[END_REF]. The initial concentration of prothrombin C II = 933.0 nM is not infered.

	tial concentration. A good trend is however observed in
	the reference and 50% cases.

Acknowledgements We are grateful to the CONACyT, Mexico scholarship and the LabEx Numev (convention ANR-10-LABX-20) for their financial support. This work was performed using HPC resources from GENCI-CINES (Grants 2016-c2015037194 and 2018-A0040307194) and with the support of the High Performance Computing Platform HPC@LR.

Conflict of Interest:

The authors declare that they have no conflict of interest.


Zihao Wang
email: zihao.wang@inria.fr

Member, IEEE Zhifei Xu

Jiayi He
email: hejiay@mst.edu

Hervé Delingette
email: herve.delingette@inria.fr

Fellow, IEEE Jun Fan
email: jfan@mst.edu

Long Short-Term Memory Neural Equalizer

Keywords: decision feedback equalizer, neuromorphic computing, deep learning, LSTM

A trainable neural equalizer based on Long Short-Term Memory (LSTM) neural network architecture is proposed in this paper to recover the channel output signal. The current widely used solution for the transmission line signal recovering is generally realized through DFE or FFE-DFE combination. The novel learning-based equalizer is suitable for highly non-linear signal restoration thanks to its recurrent design. The effectiveness of the LSTM equalizer is shown through an ADS simulation channel signal equalization task including a quantitative and qualitative comparison with an FFE-DFE combination. The LSTM neural network shows good equalization results compared to that of the FFE-DFE combination. The advantage of a trainable LSTM equalizer lies in its ability to learn its parameters in a flexible manner, to tackle complex scenario without any hardware modification. This can reduce the equalizer implantation cost for variant transmission channels and bring additional portability in practical applications.

I. INTRODUCTION

N OWADAYS, the breakthrough of massive data process- ing and data centers increases the demand for highly efficient data transmission at macro-level (telecommunication), meso-level (onboard systems), and micro-level (inside chips). With the dramatic increase of transmission speed, it is challenging to maintain the signal integrity (SI), power integrity (PI) in high-speed links design. The impedance discontinuities, field interaction, simultaneous switching noise (SSN), power delivery networks (PDN), etc, will distort the high-speed signals [START_REF] Fan | Signal integrity design for high-speed digital circuits: Progress and directions[END_REF], [START_REF] Kim | Special issue on pcb level signal integrity, power integrity, and emc[END_REF], [START_REF] Xu | Direct time-domain tan model of 3d multilayer hybrid pcb: Experimental validation[END_REF], will increase the intersymbol interference (ISI) to a point where the output of a high-speed channel may hardly be recognized by the receiver. This is why equalization techniques have been developed to reconstruct the signal at the output of the high-speed channel.

A typical Serializer-Deserializer (SerDes) system targeted for high-speed data transmission includes three modules (see Fig. 1): transmitter TX, receiver RX, and phase-locked loop (PLL) block [START_REF] He | An 8.5-12.5ghz multi-pll clock architecture with lc pll and ring pll for multi-lane multi-protocol serdes[END_REF]. The signal is transmitted from TX through the channel and received by the receiver. The PLL performs clock signal synchronization which is a naive idea for identifying the delay of the transmission channel that properly match with the delayed signal. In the RX block, the equalization techniques are implemented to reconstruct the deformed signal due to the channel loss and other complex environmental influence. There are three types of equalizers that are generally employed in the RX module. The feed-forward equalizer (FFE), the decision feedback equalizer (DFE) [START_REF] Wang | A 21-gb/s 87-mw transceiver with ffe/dfe/analog equalizer in 65-nm cmos technology[END_REF], [START_REF] Ahn | An 18-gb/s nrz transceiver with a channel-included 2-ui impulse-response filtering ffe and 1-tap dfe compensating up to 32-db loss[END_REF], [START_REF] Zheng | A 50-112-gb/s pam-4 transmitter with a fractional-spaced ffe in 65-nm cmos[END_REF], and the continuoustime linear equalizer (CTLE) [START_REF] Li | Ctle adaptation using deep learning in high-speed serdes link[END_REF]. The FFE is dedicated to compensating the loss of the channel which is assumed to be an FIR filter with a transfer function inverse to the channel transfer function. Since the compensation of that loss only is not sufficient due to the noise from the complex environment, the FFE and DFE are jointly used to have a better performance on the signal reconstruction. Since all bit signals are pulses, when the bit signal goes through the channel, there will be delays and tails. The DFE is applied to drop the tail down to 0, then the pulse tail will not affect the next bit signal, thus, the better signal reconstruction can be reached. However, due to the design of the channel and inherent non-linear representation limitations, the FFE and DFE combinations are unable to equalize original signals especially with the complex channel in very high-speed data transmission. Model based equalizers usually requires a good understanding of the transmission systems attributions and a rich expert knowledge in design but often lack of generalization ability [START_REF] Zhou | Adann: Adaptive neural network-based equalizer via online semi-supervised learning[END_REF]. Facing this challenge with the emergence of artificial intelligence technology [START_REF] Tohidi | Submodularity in action: From machine learning to signal processing applications[END_REF], several learning-based equalizers have been proposed that are more adaptable over varied channel conditions. Most of the learning-based equalizers are nonlinear equalizers. Some representative works [START_REF] Chen | Adaptive bayesian equalizer with decision feedback[END_REF], [START_REF] Chern | Semi-blind channel estimation scheme with bayesian dfe for prp-ofdm system[END_REF], [START_REF] Kiguradze | Bayesian optimization for high-speed channel equalization[END_REF] were using a Bayesian statistical model for channel equalization. Neural network-based approaches are used by some prior works [START_REF] Abdulkader | Neural networks-based turbo equalization of a satellite communication channel[END_REF], [START_REF] Rajbhandari | Neural network-based joint spatial and temporal equalization for mimovlc system[END_REF], [START_REF] Zhang | 137 gb/s pam-4 transmissions at 850 nm over 40 cm optical backplane with 25 g devices with improved neural network-based equalization[END_REF], [START_REF] Zhou | Adann: Adaptive neural network-based equalizer via online semi-supervised learning[END_REF] for channels equalization. Some uncommonly used low-speed equalization technologies such as Support Vector Machine based equalizer [START_REF] Taub | Svm based method for multi-equalizer optimization[END_REF], [START_REF] Bi | A powerful equalizer based on modified svm classifier without nonlinear kernel enabled 100-gb/s ng-epon system with 10-g class[END_REF], Fuzzy based networks equalizer [START_REF] Mayer | Blind fuzzy adaptation step control for a concurrent neural network equalizer[END_REF] were more difficult for hardware implementation and are not elaborated in detail here.

Recent work of Chu et al. [START_REF] Chu | Equalization with neural network circuitry for high-speed signal link[END_REF] introduced a neural network circuitry equalizer (NNE) for high-speed signal link construction. An artificial neural network (ANN) is integrated into the continuous-time linear EQ (CTLE) as the replacement of the DFE block (block (b) in Fig. 1). Their quantitative experiments have shown that the ANN-based equalizer is more robust than the DFE approach for degraded signal restoring accuracy. In terms of the transmission speed improvement, the NNE improved about 30% in data rate vis-à-vis the traditional equalization method. A similar work of Rajbhandari et al. employed ANN for multiple-input-multiple-output visible light communication (MIMO-VLC) system signal equalization. The MIMO-VLC equalization problem is usually more challenging than the PCB level signal equalization as the channel of the MIMO-VLC system can be an open environment which may be easily violent by the external optical noisy. The experimental study has shown that the ANN equalizer outperformed the temporal-spatial cross-talks compensation [START_REF] Rajbhandari | Neural network-based joint spatial and temporal equalization for mimovlc system[END_REF] and possible for the hardware implementation.

The conventional ANN-based equalizer is not capable of capturing the time series information as the input signals are processed without considering the temporal relationship [START_REF] Dai | A comparison of mlp, rnn and esn in determining harmonic contributions from nonlinear loads[END_REF]. However, the signal equalization task is typically a temporal processing problem that is suitable for employing our LSTM based equalizer. Several prior works have shown that the high potential of neural network-based equalization for the future high-speed, complexity electromagnetic channels and longdistance transmission systems design. We see that most of the neuromorphic computing systems used for equalization are still limited to the basic perceptions or multi-layer perceptron neural systems.

In this paper, we propose a novel approach based on the long short-term memory neural network and deep learning for channel signal equalization. In many circumstances, neural network computing is software implemented and run on GPGPU (General-purpose processing on graphics processing units) or TPU (Tensor Processing Unit) acceleration which limit its application for highly integrated circuit systems. Yet, the proposed approach can be implemented through programmable gate array or ASIC or general neuromorphic computing instruments (related literature: [START_REF] Rajendran | Low-power neuromorphic hardware for signal processing applications: A review of architectural and system-level design approaches[END_REF], [START_REF] Lecun | 1.1 deep learning hardware: Past, present, and future[END_REF], [START_REF] Balaji | Mapping spiking neural networks to neuromorphic hardware[END_REF]) which realize the embedding for the integrated system. To the best of our knowledge, our proposed method is the first application of a physical realizable temporal memory-based neural equalizer. Furthermore, the proposed approach is adaptable to variable frequency signal equalization thanks to the ability of neural networks for performing powerful approximations.

The paper is organized as follows: In section II we first describe the memory gates mechanism of the LSTM network and the construction procedures for LSTM based equalizer architecture design. Section III describes the implementation of the proposed method both simulated and physical realization. The conclusion section V summarizes the full paper and look into the future of applying deep learning approaches for equalizer design.

II. LONG SHORT-TERM MEMORY NEURAL EQUALIZER

A. Long Short-Term memory neural network

The long Short-Term memory neural network (LSTM) [START_REF] Hochreiter | Long short-term memory[END_REF] is a type of recurrent neural networks (RNN) which is specially designed for processing time series or sequential data. Differing from fully connected neural networks, RNNs can learn temporal information by performing recurrent backpropagation [START_REF] Rumelhart | Learning Representations by Back-propagating Errors[END_REF]. Fig. 2 shows a RNN for signal series regularization. At time t, the RNN is fed with a noisy signal x t = r t and the memory information R t-1 passed from the last time t-1. The neural network h t will then generate as output the processed signal O t = s t at time t and memory information R t . Yet, RNN is limited by the transmission of past information since the memory R t-1 might be not sufficient for processing the input of current time t. The LSTM design tackles this problem by introducing the long short-term memory gates mechanism [START_REF] Hochreiter | Long short-term memory[END_REF]. It has been shown in practice that the LSTM is better than RNN for long-distance memory processing, and LSTM has been successfully applied to many areas related to time series signals such as Google neural machine translation system, Microsoft voice recognition system, ECG beats signal classification, and Apple's Siri assistant systems, etc [START_REF] Wu | Google's neural machine translation system: Bridging the gap between human and machine translation[END_REF], [START_REF] Li | Lstm time and frequency recurrence for automatic speech recognition[END_REF], [START_REF] Übeyli | Combining recurrent neural networks with eigenvector methods for classification of ecg beats[END_REF].

x t-1 h t-1 o t-1 U W x t h t o t U W x t+1 h t+1 o t+1 U W R R R
As shown in Fig. 3, the LSTM cell consists of a group of neurons and several information transformation flows. The cell state variable c t-1 carries the memory information of the previous cell whereas the hidden variable h t-1 is the cell output information of the previous cell. Within the LSTM cell, two activation functions are used to perform non-linear transformations of the information flows: • The Sigmoid function:

σ(x) = 1 1+exp(-x) ∈]0, 1[ com- putes σ(h t-1
) which controls the input gate I t , and the forget gate F t by a Hadamard product (or element-wise multiplication noted ) with the last cell state c t-1 . The activation level of the forget gate is correlated with the information passing from the last cell state.

• The hyperbolic tangent function :

tanh(x) = exp(2x)-1 exp(2x)+1 ∈] -1,
f t = σ(w f x t + w rf h t-1 + b f ) i t = σ(w i x t + w ri h t-1 + b i ) o t = σ(w o x t + w ro h t-1 + b o ) cs t = tanh (w c x t + w rc h t-1 + b c ) (1)
where f , i, o are the forget , input , and output gates while cs is the cell input activation vector. The update rules for the memory (c t ) and output (h t ) are given by :

c t = f t c t-1 + i t cs t h t = o t tanh(c t ) (2)
Thus, the current cell state is a combination of the previous cell state and the cell input activation controlled by the forget and input gate variables. The hidden state is a function of the cell state weighted by the output gate variable.

In addition, the trainable parameters for the LSTM networks are:

• The neural networks weight parameters:

W = {{w n f , w n i , w n c , w n o ]}, n ∈ N } (3)
• The bias parameters :

b = {{b n f , b n i , b n c , b n o ]}, n ∈ N } (4)
• The recurrent weights:

W r = {{w n rf , w n ri , w n rc , w n ro ]}, n ∈ N } (5)
Furthermore, an LSTM cell can be easily implemented in a digital circuit with adjustable neural activation parameters and the LSTM circuit is amenable to analog computing in which the neurons are directly modeled by electronic signals [START_REF] Han | Era-lstm: An efficient reram-based architecture for long short-term memory[END_REF], [START_REF] Zhao | Long short-term memory network design for analog computing[END_REF], [32]. The analog neural network implementation can overcome the frequency limitation of equalizers where the main issue is the frequency of the converter of the digital and analog signal (DAC).

B. Long short-term memory neural equalizer

, c t-1 } f t ← σ(w f r t + w rf h t-1 + b f ) ; // forget gate i t ← σ(w i r t + w ri h t-1 + b i ) ; // input gate cs t ← tanh (w c r t + w rc h t-1 + b c ) ; // cell gate o t ← σ(w o r t + w ro h t-1 + b o ) ; // output gate // Compute the current cell outputs: c t , h t c t ← f t ⊗ c t-1 + i t ⊗ cs t ; // cell activation h t ← ø t ⊗ tanh(c t ); // hidden gate h, c ← h t , c t ;
// Refresh states cache s * ← F(h); // Fully connected layer ŝ ← FIR(s * ); // Post processing filtering end

The long short-term memory neural equalizer is an intrinsicfeedback equalizer system where the spoiled signal series r shown in Alg: 1 from the transmission line are directly processed by the LSTM neural block while the feedback information is achieved by the long-short memories states h, c. As Fig. 4 shows, the LSTM is integrated as a terminal The output hidden state h t at time t is decoded by a fully connected neuron F(h t ) to get the equalized signal output. This fully connected cell is parameterized by w f ull and b f ull . The output of the full connected cell is given by:

F(h t ) = σ(w f ull ḣt + b f ull ) (6)
Finally, the signal is processed by a FIR filter to remove the noise and stabilize the waveform. The LSTM network aims to minimize the difference between the equalized bits ŝ and the bits signal s emitted from the chip (a). We can train the LSTM by minimizing the Mean-Square-Error (MSE) of the two items:

M SE = N i-1 (s i -s i ) 2 N (7)
The LSTM network parameters can be trained through the backpropagation algorithm [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF] to estimate the gradient of the 8) represents a sequence of signals which terminates at time t (see Fig. 5), and the same notation is used for the validation dataset D K valid : {S valid |R valid }. The training dataset can be either measured from physical transmission platforms or can be collected from a simulation of the transmission channels when the S-parameters of the transmission line is known. In the training phase, the D K train is fed to the LSTM network by batch to minimize the MSE loss function by stochastic gradient descent based algorithms (i.e.: Adam in this work) [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF].

R train =        {R t-n . . . R t } 1 {R t-n-1 . . . R t-1 } 2 . . . {R 0 . . . R n } K        ; S train =        {S t } 1 {S t-1 } 2 ... {S n } K        (8) R t-n ...R t in equation (
b) Backpropagation: The parameters of the LSTM equalizer are determined through a batch-based backpropagation algorithm [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF] which is performed in two steps: (I) forward propagation as described in Alg: 1. (II) backward propagation where the neural network parameters are updated based on the gradients' information of error feedback. The details of backpropagation are described in Alg: 2

III. EXPERIMENTS AND RESULTS

A. Experimental architecture

The proposed LSTM equalizer proposed in this paper was implemented and verified through Simulink modeling for a

B. Experimental result for a high-speed PCB channel

A PCB channel, containing transmission line, via, connector, etc, has been used for this experiment. The S parameter of the channel characterized up to 50 GHz is shown in Fig. 8 with both insertion and return losses. The insertion loss exposes very high loss in high frequencies. The transmitted signal is severely attenuated and therefore hard to recognize by the receiver. Therefore, a FFE-DFE or our proposed method is needed to recover the signal.

A simulation is performed in ADS with the S parameters of the fabricated PCB channel. The schematic is shown in Fig. 9. There are two channels within this S-parameter file, one of the channels is employed in this simulation for simplicity. A 50 Gbps data rate is chosen for the simulation. The high/low levels are defined as 1 V /0 V , respectively. The rise and fall times are equal to 4 psec. // Initialize states: h = 0, c = 0; // while valid < δ converge do for i = 1; i < K; i

+ + do R i ← {R t-n ...R t } i ; // Get ith input S i ← {S t } i ; // Get ith training target // Forward propagation with Algo: 1 R ← LST M (R i );
// // Get the loss between output and target based on loss function:7

= loss(R , S i ); // ∆ W ← α × ∂ ∂W ; // Jacobian W ∆ W R ← α × ∂ ∂W R ; // Jacobian WR ∆ b ← α × ∂ ∂b ; // Jacobian b W ← ∆ W + W ; // Update W W R ← ∆ W R + W R; // Update WR b ← ∆ b + b; // Update b if 0 ≡ i mod V then valid ← 0 ; // for j = 1; j < K ; j + + do // Get validation dataset R j ← D K valid : {R t-n ...R t } j ; // S j ← D K valid : {S t } j ; // // Forward propagation with Algo: 1 R ← LST M (R j );
// valid ← valid + loss(R , S j ); // end valid = valid K ; // break; // else continue; end end end values applied in this case are shown in Table I. The active FFE is applied to recover the loss due the channel physical parameters.

Yet, with the FFE-DFE equalizer processing, the eye is opened slightly as shown in Fig. 11 with eye height is 0.18 V, eye width is 5.1 ps and jitter is 13.8 ps.

Fig. 12 shows the eye diagram when using the LSTM equalizer. The eye characteristics comparison between the active FFE-DFE and LSTME is shown in Table . II. We see that the eye width, height and jitter are much improved compared to the active FFE-DFE approach.

2) Waveform comparison: The waveform comparison is shown in Fig. 13. The red curve is the output signal without any equalization whereas the blue signal in the top (resp. bottom) figure is the output signal after the LSTM (resp. active FFE/DFE) equalization. Even with active FFE-DEF equalization, the loss and shape of the waveform are significantly degraded compare to the original signal. However, the case

C. Experimental results for a mismatched channel

In the world of signal integrity, due to the limitation of PCB manufacturing, impedance mismatching is a general problem for high-speed transmission system designs. The overshoot and undershoot problems will appear due to the impedance mismatching. And the defect on the channel will also distort the output signals. Therefore, in this experiment, our proposed equalizer is used for a badly designed channel with mismatched loads to exhibit the strength of our method. The channel characteristic is shown in Fig. 14, the channel exhibit the high reflection which is mainly due to the impedance mismatching. The eye diagram for this badly designed channel is first obtained through ADS simulation with optimized FFE and DFE settings. A 20 Gbps signal with 10 ps rise/fall time is applied to this badly designed channel. A limitation of LSTM equalizer lies in the relative increase of computational complexity compared to the traditional FFE-DFE combination. This increase is related to highly non-linear representation capability of deep neural networks. Empirically, we noted that in section III for the equalization of a fabricated PCB channel with a fixed frequency 50 Gbps input signal, an LSTM equalizer with 20 hidden cells and 15 delays with the resolution of 5 ps is sufficient to get a much improved equalized output. In this case, the computational complexity is acceptable for practical use. However, the complexity of a deep learning LSTM equalizer may increase significantly with the demand for variable-frequency equalization. Therefore, a more thorough study of the required neural model complexity with respect to variable equalization frequency, and different transmission channels model is necessary.

The proposed LSTM equalizer is based on supervised learning with a consistent paired training dataset to learn the non-linear mapping. This requires to collect specific simulated data and measurements, although it is also necessary for FFE-DFE based channel equalization. With the recent advances of semi-supervised/unsupervised learning [START_REF] Zhai | S4l: Self-supervised semi-supervised learning[END_REF], [START_REF] Jing | Self-supervised visual feature learning with deep neural networks: A survey[END_REF], there is a possibility to reduce the efforts for collecting the training dataset by applying those self-adaptive learning approaches to signal equalizer design.

V. CONCLUSION

In this paper, we have proposed a novel recurrent feedback equalizer based on the LSTM neural network. The LSTM network is trained offline with a collected signal sequence dataset. We have shown that our approach outperforms commonly used FFE-DFE signal equalization methods with meaningful metrics and results. Furthermore, the LSTM equalizer is amenable to hardware/software implementation with flexible portability for different channels signals equalization. Future research will address the issue of variable frequency signal equalization and weakly-supervised learning approaches for the next-generation AI-based equalizer design. Moreover, the compatibility of the neural learning-based approaches for high-frequency, low-cost, low-power, and low-latency signal equalization needs to be further studied.

Fig. 1 :

 1 Fig. 1: Backboard SerDes architecture. Signal emission from chips is transferred through a transmission line which leads to loss and distortion.

Fig. 4 :

 4 Fig. 4: The signal link of the proposed LSTM neural network based channel equalizer. (a) Chip generates logic bits. (b) Arbitrary signal channel. (c) N bits signal delay. (d) LSTM neural network signal interpreter. (e) Signal receiver.

Fig. 5 : 15 Fig. 6 :

 5156 Fig. 5: LSTM Equalizer diagram shows input spoiled signals R t and equalized signals S t . The cell states and hidden states are feedback signals between different time-frame.

  MSE functional on a batch of data. C. LSTM training a) Dataset: Both a training and a validation dataset need to be collected for training the LSTM equalizer. The neural network is trained through a batch B t of training data and corresponding target signal. The validation dataset is used to control the training process in order to avoid overfitting the training dataset. We write the training set as D K train : {S train |R train } in Alg: 2, where

Fig. 7 :

 7 Fig. 7: The Simulink simulation architecture diagram of the proposed LSTME applied for a simulated channel signal equalization.

Fig. 8 :

 8 Fig. 8: Return loss and insertion loss of the simulated channel.

Fig. 9 :

 9 Fig. 9: Schematic of the ADS channel simulation. 1) Eyediagram comparison: To assess the performance of the proposed approach, we use eye diagram analysis which provides intuitive evaluation for the performance of digital signal equalization. As shown in Fig. 10, the eye is not open at all in the eyediagram of the output data which means the output signals will have high bit-error rate (BER) without the FFE-DFE equalization. The FFE and DFE cursors

Fig. 10 :

 10 Fig. 10: Eye diagram of ADS channel signal simulation (raw signal from channels).

Fig. 11 :

 11 Fig. 11: Eye diagram of ADS channel signal equalized via the conventional active FFE-DFE system.

Fig. 12 :

 12 Fig. 12: Eye diagrams of ADS channel signal (50 GHz) equalized via the proposed LSTM equalizer system.

Fig. 13 :

 13 Fig. 13: Comparison between LSTME (top) and FFE-DFE (bottom) equalized signal waveform with input and output signal of the transmission channel shown.

Fig. 14 :

 14 Fig. 14: Return loss and insertion loss of the bad designed channel

Fig. 15 :

 15 Fig. 15: ADS simulated eye diagram with bad designed channel

  . The blue curve shown in the lower figure of Fig. 17 indicates the output waveform of the active FFE-DFE. The waveform is significantly degraded with overshoot and undershoot. However, for the output waveform of the LSTME shown in the top figure of Fig.

17

 17

  exposes perfect match to the original signal (red curve). The efficiency of our proposed method is proved again.IV. DISCUSSIONThe output waveforms with different equalization techniques are shown in Fig.13and Fig.17

Fig. 16 :

 16 Fig. 16: Eye diagrams of badly designed channel signal (20 GHz) equalized via the proposed LSTM equalizer system.

Fig. 17 :

 17 Fig. 17: Comparison of the waveform quality of the LSTM equalizer (Top) and FFE-DFE (Bottom) for a badly designed channel of 20 GHz signal.

  1[ is used for memory selection by controlling the activation level of the current cell input state x t . The computation of the intermediate variables shown in Fig: 3 are given by:

Algorithm 1 :

 1 Signal forward propagation at time t.

		// Get last cell
	states
	r t ← r n mem ;	// Get current signals
	// Compute the three cell gates and
	activation based on the input
	signal vector r t and last states
	{h t-1

Inputs: N bit Signal delayed: r n mem Output: Equalized signal: s Restore neural parameters from RAM: W, W R, b; // LSTM Initialize states cache: h = 0, c = 0; // while Clock Refresh do t ← clock signal h t-1 , c t-1 ← h, c ;

TABLE I :

 I DFE and FFE parameters

	FFE	DFE
	PreCursor	PostCursor	Taps
	-2.337340	4.038660	0.322812
			-0.017401
		-2.185680	0.048581
	0.782150	0.534350	-0.065590
			0.039204
		-0.121820	-0.021085

TABLE II :

 II Eyediagram characteristic comparison

		Jitter (ps) Eye width (ps) Eye height (V)
	Active FFE-DFE	13.8	5.1	0.18
	LSTME	9.8	12.8	0.87

Manuscript received xxx. This work has been partially supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. And in part by the US National Science Foundation under Grant No. IIP-1916535. Zihao Wang and Hervé Delingette are with the French National Institute for Research in Computer Science and Control (INRIA) and Université Côte d'Azur (UCA), and Hervé Delingette is with French Interdisciplinary Institutes of Artificial Intelligence (3IA Côte d'Azur). 06902 Valbonne,