Rodrigo Méndez Rojano

| Mansur Zhussupbekov

James F Antaki

| Didier Lucor

Uncertainty quantification of a thrombosis model considering the clotting assay PFA-100®

Keywords: Thrombosis, PFA-100®, uncertainty quantification, polynomial chaos expansion, sensitivity analysis

Mathematical models of thrombosis are currently used to study clinical scenarios of pathological thrombus formation. As these models become more complex to predict thrombus formation dynamics high computational cost must be alleviated and inherent uncertainties must be assessed. Evaluating model uncertainties allows to increase the confidence in model predictions and identify avenues of improvement for both thrombosis modeling and anti-platelet therapies. In this work, an uncertainty quantification analysis of a multi-constituent thrombosis model is performed considering a common assay for platelet function (PFA-100®). The analysis is facilitated thanks to time-evolving polynomial chaos expansions used as a parametric surrogate for the full thrombosis model considering two quantities of interest; namely, thrombus volume and occlusion percentage. The surrogate is thoroughly validated and provides a straightforward access to a global sensitivity analysis via computation of Sobol' coefficients. Six out of fifteen parameters linked to thrombus consitution, vWF activity, and platelet adhesion dynamics were found to be most influential in the simulation variability considering only individual effects; while parameter interactions are highlighted when considering the total Sobol' indices. The influential parameters are related to thrombus constitution, vWF activity and platelet to platelet adhesion dynamics. The surrogate model allowed to predict realistic PFA-100® closure times of 300,000 virtual cases that followed the trends observed in clinical data.

The current methodology could be used including common anti-platelet therapies to identify scenarios that preserve the hematological balance.

INTRODUCTION

Thrombosis, which is defined as excessive formation of blood clot or thrombus, is a common pathology in several cardiovascular diseases [START_REF] Nagareddy | Inflammation and thrombosis in cardiovascular disease[END_REF] , and blood-wetted medical devices [START_REF] Dal Sasso | The Biocompatibility Challenges in the Total Artificial Heart Evolution[END_REF][START_REF] Jaffer | The blood compatibility challenge. Part 1: Blood-contacting medical devices: The scope of the problem[END_REF][START_REF] Labarrere | Thrombogenic and Inflammatory Reactions to Biomaterials in Medical Devices[END_REF] . Thrombus formation is characterized by an intertwined process of platelet activity and coagulation. In this hemostatic process, platelet activation and aggregation leads to formation of a platelet plug that is mechanically stabilized by a fibrin net formed by coagulation reactions [START_REF] Colman | Hemostasis and Thrombosis[END_REF] . In hemodynamic conditions that promote high shear stresses, such as stenotic vessels or prosthetic heart valves or blood pumps, von Willebrand Factor (vWF) plays a central role in thrombosis. vWF is a blood glycoprotein that has a collapsed globular conformation in its natural state and unfolds in response to extensional flow and high-shear flow conditions. Once unfolded, vWF plays an important role in the platelet adhesion and aggregation process enabling clot formation.

Multi-scale computational models of thrombosis have been developed in the past decades to understand and predict the thrombus formation dynamics in academic and clinical configurations [START_REF] Leiderman | Grow with the flow: A spatial-temporal model of platelet deposition and blood coagulation under flow[END_REF][START_REF] Menichini | Predicting false lumen thrombosis in patient-specific models of aortic dissection[END_REF][START_REF] Yazdani | A General Shear-Dependent Model for Thrombus Formation[END_REF][START_REF] Xu | A multiscale model of thrombus development[END_REF][START_REF] Zhiliang | Fibrin Networks Regulate Protein Transport during Thrombus Development[END_REF][START_REF] Xu | Model predictions of deformation, embolization and permeability of partially obstructive blood clots under variable shear flow[END_REF][START_REF] Yazdani | Data-driven Modeling of Hemodynamics and its Role on Thrombus Size and Shape in Aortic Dissections[END_REF] ). These models span several degrees of complexity, from simple flow characteristic indices to complicated, coupled biochemical interactions. Mechanistic models, for example, may contain equations that describe shear dependent platelet activity including adhesion, aggregation, and activation [START_REF] Sorensen | Computational simulation of platelet deposition and activation: I. Model development and properties[END_REF][START_REF] Storti | dFN. A continuum model for platelet plug formation and growth[END_REF][START_REF] Danes | A density-dependent FEM-FCT algorithm with application to modeling platelet aggregation[END_REF] ; hemodynamics that regulate the transport of biochemical species [START_REF] Storti | A continuum model for platelet plug formation, growth and deformation[END_REF][START_REF] Fogelson | Fluid mechanics of blood clot formation[END_REF] ; and/or coagulation reactions that lead to the formation of fibrin [START_REF] Anand | A model for the formation and lysis of blood clots[END_REF][START_REF] Yesudasan | Recent advances in computational modeling of fibrin clot formation: A review[END_REF][START_REF] Hansen | Automated reduction of blood coagulation models[END_REF] . Increasing modeling complexity improves the fidelity and versatility of the model, but the computational cost of the simulation escalates. Also, increasing detail introduces additional uncertainties that can have a great impact on its accuracy [START_REF] Belyaev | Modeling thrombosis in silico: Frontiers, challenges, unresolved problems and milestones[END_REF] . These uncertainties are related to the fundamental structure of the model , numerical discretization, and/or model parameters such as diffusion coefficients, concentrations of biochemical species, fluid viscosity, etc. Previous investigators have introduced several uncertainty quantification (UQ) strategies to evaluate the effect of such uncertainties on model performance in the context of blood flows [START_REF] Sankaran | Uncertainty quantification in coronary blood flow simulations: Impact of geometry, boundary conditions and blood viscosity[END_REF][START_REF] Brault | Uncertainty quantification of inflow boundary condition and proximal arterial stiffness-coupled effect on pulse wave propagation in a vascular network[END_REF][START_REF]Uncertainty Quantification of a Multiscale Model for In-Stent Restenosis[END_REF][START_REF] Steinman | Editorial: Special Issue on Verification, Validation, and Uncertainty Quantification of Cardiovascular Models: Towards Effective VVUQ for Translating Cardiovascular Modelling to Clinical Utility[END_REF] . UQ analyses in thrombosis models have been performed to identify the sensitivity of input parameters [START_REF] Danforth | The impact of uncertainty in a blood coagulation model[END_REF][START_REF] Rojano | Kinetics of the coagulation cascade including the contact activation system: sensitivity analysis and model reduction[END_REF][START_REF] Link | A mathematical model of coagulation under flow identifies factor V as a modifier of thrombin generation in hemophilia A[END_REF][START_REF] Melito | Sensitivity Analysis of a Phenomenological Thrombosis Model and Growth Rate Characterisation[END_REF] . In most sensitivity studies, the global Morris method [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF] or the Sobol' indices method [START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF] have been used. The Morris method is used as a screening tool to identify important parameters in models with a large number of input variables. The Sobol' method produces indices based on variance that quantifies the output uncertainty related to each input parameter. For example, Melito et al. [START_REF] Melito | Sensitivity Analysis of a Phenomenological Thrombosis Model and Growth Rate Characterisation[END_REF] performed a variance based sensitivity analysis on the thrombosis model of Menichini et al. [START_REF] Menichini | Predicting false lumen thrombosis in patient-specific models of aortic dissection[END_REF] considering thrombosis in a backward facing step. Their results suggest that only four of the nine parameters included in their thrombosis model had significant impact in thrombus formation implying that the model could be made sparser to optimize computational resources. Link et al. [START_REF] Link | Computationally Driven Discovery in Coagulation[END_REF] performed a sensitivity analysis to identify important coagulation factors in the thrombin generation profile among hemophilia A patients. Their work allowed the identification of coagulation factor V as an important modifier of thrombin formation among patients with hemophilia A, paving the way to pharmacological interventions to treat these patients. In addition, as demonstrated by Link et al. [START_REF] Link | Computationally Driven Discovery in Coagulation[END_REF] the synthetic databases built to perform sensitivity analysis can serve as hypothesis generation tool for in-vitro experiments or, in the long run, patient treatment.

In this work, a forward uncertainty propagation study is conducted using a well-established but computationally costly multiconstituent thrombosis model introduced by Wu et al. [START_REF] Wu | Multi-Constituent Simulation of Thrombus Deposition[END_REF] . A widely used platelet function assay PFA-100® (Siemens, Erlangen, Germany) was selected as the study case. To understand the impact of uncertain input parameters a variance based global sensitivity analysis was conducted monitoring Sobol' sensitivity coefficients taking advantage of a polynomial chaos approximation as a surrogate of the baseline thrombosis model. The surrogate model is also used in this study to assess its predictive ability of the PFA-100® closure time distribution in three scenarios: normal platelet count, thrombocytopenia, and vWF deficiency. The predicted closure times distribution align well with clinical data.

METHODS

Baseline Thrombosis Model

The thrombosis model used in the current work is a variation of the approach introduced by Wu et al. [START_REF] Wu | Multi-Constituent Simulation of Thrombus Deposition[END_REF] with the addition of vWF activity. Figure 1 A shows how platelets, the central player of the model, deposit to artificial surfaces. Figure 1 B illustrates vWF unfolding due to high shear rates increasing platelet deposition. Thrombus growth is quantified through platelet deposition over artificial surfaces and subsequent platelet aggregation. Several biochemical species are considered in the model: 1) AP Activated platelets, 2) RP Resting Platelets, 3) AP 𝑑 Deposited activated platelets, 4) RP 𝑑 Deposited Resting Platelets 5) ADP Adenosine Diphosphate, 6) TxA 2 Thromboxane A2, 7)TB Thrombin , 8) PT Prothrombin, 9) AT Anti-Thrombin, 10) vWF c vWF Coiled, and 11) vWF s vWF stretched.

The spatial and temporal dynamics of biochemical species are quantified through convection-diffusion-reaction equations of the form:

𝜕𝑐 𝑖 𝜕𝑡 = ∇ ⋅ (𝐷 𝑖 ∇𝑐 𝑖 ) -𝒗 𝑓 ⋅ ∇𝑐 𝑖 + 𝑟 𝑖 , (1)
where: 𝑐 𝑖 is the concentration of species 𝑖, 𝐷 𝑖 is the diffusion coefficient, 𝒗 𝑓 is the velocity vector field, and 𝑟 𝑖 is the reaction source term that accounts for biochemical interactions such as platelet activation by thrombin or TxA 2 .

vWF activity

To account for the role of vWF in high shear stress thrombosis the original model of Wu et al. [START_REF] Wu | Multi-Constituent Simulation of Thrombus Deposition[END_REF] was amended to include the enhanced platelet aggregation and adhesion driven by vWF activity. The right-hand side of Figure 1 shows the two interconvertible states of vWF in the model: collapsed vWF c and stretched vWF s . In all simulations, only vWF c is introduced at the inlet of the domain and vWF s is produced by means of vWF c unfolding. The presence of stretched vWF s polymers has a local thrombogenic effect via two mechanisms: it amplifies the deposition rate of free-flowing platelets by scaling it by 𝛼 𝑘 (Eq (2)) and increases the resistance of deposited platelets to shear cleaning by scaling the shear-cleaning parameters, 𝜏 𝑒𝑚𝑏 and 𝜏 𝑒𝑚𝑏,𝑏 , by 𝛼 𝜏 according to Eq (3)

𝛼 𝑘 = { 1, if vWF s ≤ vWF crit s vWF s vWF crit s , if vWF s > vWF crit s (2) 𝛼 𝜏 = { 1, if vWF s ≤ vWF crit s 𝐷 𝑒 ( ̇𝛾) 𝐷 𝑙 𝑒 , if vWF s > vWF crit s , (3)
where vWF crit s is set as 5% of the total vWF concentration; 𝐷 𝑒 is the Morse potential's well depth proposed by Yazdani et al. 8 that correlates the platelet adhesion force to the local shear rate, ̇𝛾, according to Eq (4);

𝐷 𝑒 ( ̇𝛾) = 𝐷 ℎ 𝑒 [ tanh ( ̇𝛾 -̇𝛾𝑣𝑊𝐹 1000 ) + 𝐷 𝑙 𝑒 𝐷 ℎ 𝑒 + 1 ] , ( 4
)
𝐷 𝑙 𝑒 and 𝐷 ℎ 𝑒 determine the adhesive forces at low and high shear rates, respectively, and ̇𝛾𝑣𝑊𝐹 = 5500 s -1 is the critical shear rate for vWF unfolding and marks the transition from low-shear to high-shear regime in this expression. The value for 𝐷 𝑙 𝑒 , 𝐷 ℎ 𝑒 , vWF crit s , and ̇𝛾𝑣𝑊𝐹 , were obtained from Yazdani et al. [START_REF] Yazdani | A General Shear-Dependent Model for Thrombus Formation[END_REF] . The scalar shear rate ̇𝛾 is computed following OpenFoam source code 34 :

̇𝛾 = √ 2𝑇 𝑟(𝑠𝑦𝑚𝑚(∇𝒖) ⋅ 𝑠𝑦𝑚𝑚(∇𝒖)) (5)
vWF c unfolding occurs via two mechanisms: periodic tumbling in simple shear and strong unfolding in flows with extensional kinematics. The flow classification is based on polymer unfolding criteria by Babcock et al. [START_REF] Babcock | Visualization of molecular fluctuations near the critical point of the coil-stretch transition in polymer elongation[END_REF] where they define the flow-type 𝜆 as:

𝜆 = ‖𝑺‖ -‖𝜴‖ ‖𝑺‖ + ‖𝜴‖ (6)
where ‖𝑺‖ and ‖𝜴‖ are the magnitudes of the symmetric and anti-symmetric parts of the velocity gradient tensor, respectively. Thus, in the limiting cases of the flow type parameter, 𝜆 = 1 corresponds to purely extensional flow, 𝜆 = 0 to simple shear flow, and 𝜆 = -1 to solid-body rotation. vWF will experience strong unfolding in flows with 𝜆 ≥ 𝜆 𝑐𝑟𝑖𝑡 , simple-shear tumbling in -𝜆 𝑐𝑟𝑖𝑡 ≤ 𝜆 < 𝜆 𝑐𝑟𝑖𝑡 , and remain collapsed if 𝜆 < -𝜆 𝑐𝑟𝑖𝑡 .

In simple shear and near-shear flows, -𝜆 𝑐𝑟𝑖𝑡 ≤ 𝜆 < 𝜆 𝑐𝑟𝑖𝑡 , the unfolding rate, 𝑘 𝑐-𝑠 (collapsed-to-stretched conversion), follows a function proposed by Lippok et al. [START_REF] Lippok | Shear-Induced Unfolding and Enzymatic Cleavage of Full-Length VWF Multimers[END_REF] for shear-dependent cleavage rate of vWF. Since this function describes the availability of vWF monomers for enzymatic cleavage, the unfolding rate is assumed to follow the same shape, given by Eq. ( 7)

𝑘 𝑠ℎ𝑒𝑎𝑟 𝑐-𝑠 ( ̇𝛾) = 𝑘 ′ 1 + exp ( - ̇𝛾-̇𝛾𝑣𝑊𝐹 Δ ̇𝛾 ) (7
)

𝑘 𝑠ℎ𝑒𝑎𝑟 𝑠-𝑐 = 𝑘 ′ (8)
where ̇𝛾𝑣𝑊𝐹 = 5500 s -1 is the critical (half-maximum) shear rate for vWF unfolding and Δ ̇𝛾 = 1271 s -1 denotes the width of the transition. 𝑘 ′ is the nominal state conversion rate, given by 𝑘 ′ = 1∕𝑡 𝑣𝑊 𝐹 , where 𝑡 𝑣𝑊 𝐹 = 50 ms is the vWF unfolding time. This timescale for vWF unfolding was reported in experiments by Fu et al. [START_REF] Fu | Flow-induced elongation of von Willebrand factor precedes tension-dependent activation[END_REF] and Brownian dynamics simulations by Dong et al. [START_REF] Dong | A mechano-reactive coarse-grained model of the blood-clotting agent von Willebrand factor[END_REF][START_REF] Kania | Predicting pathological von Willebrand factor unraveling in elongational flow[END_REF] To reflect the tumbling behavior of the vWF in this regime, the stretched-to-collapsed conversion rate, 𝑘 𝑠-𝑐 is set to the nominal value 𝑘 ′ (Eq. ( 8)).

In extensionally-dominated flows with 𝜆 ≥ 𝜆 𝑐𝑟𝑖𝑡 , strong unfolding occurs if the modified Weissenberg number, 𝑊 𝑖 𝑒𝑓𝑓 , given by Eq. ( 9), exceeds the critical value 𝑊 𝑖 𝑒𝑓𝑓,𝑐𝑟𝑖𝑡 . Then, the unfolding rate 𝑘 𝑠𝑡𝑒𝑎𝑑𝑦 𝑐-𝑠 scales with 𝑊 𝑖 𝑒𝑓𝑓 according to Eq. (10)

) 11
where 𝜏 𝑟𝑒𝑙 is the polymer relaxation time used to non-dimensionalize the expression. Here, we adopted the unfolding threshold and the hysteresis value reported by Sing and Alexander-Katz. [START_REF] Sing | Dynamics of collapsed polymers under the simultaneous influence of elongational and shear flows[END_REF] Assuming the monomer diffusion time of 1.02 × 10 -3 s for 𝜏 𝑟𝑒𝑙 in Eq. ( 9), 𝑊 𝑖 𝑒𝑓𝑓 ,𝑐𝑟𝑖𝑡 = 0.316 and 𝑊 𝑖 𝑒𝑓𝑓 ,ℎ𝑦𝑠𝑡 = 0.053.

If 𝜆 < -𝜆 𝑐𝑟𝑖𝑡 , vWF c remains collapsed and vWF s reverts to a globular state, so 𝑘 𝑐-𝑠 = 0 and 𝑘 𝑠-𝑐 = 𝑘 ′ . The reaction term for vWF c and vWF s is then obtained by multiplying the reaction rates previously described by the concentration of this species following the law of mass action. A detailed validation of the vWF-thrombosis model is provided in Zhussupbekov et al. [START_REF] Zhussupbekov | A Continuum Model for the Unfolding of von Willebrand Factor[END_REF] for multiple academic simple flow configurations and complex in-vitro clotting scenarios.

Flow dynamics

The pressure and velocity fields 𝑝 and 𝒗 𝑓 are obtained solving the Navier-Stokes equations considering an incompressible Newtonian fluid:

𝜕𝜌 𝑓 𝜕𝑡 + ∇ ⋅ (𝜌 𝑓 v 𝑓 ) = 0 ( 12
)
𝜌 𝑓 𝐷v 𝑓 𝐷𝑡 = ∇ ⋅ T 𝑓 + 𝜌 𝑓 b 𝑓 + 𝐶 2 𝑓 (𝜙)(v 𝑓 -v 𝑇 ) (13)
where T 𝑓 is the stress tensor of the fluid described as:

T 𝑓 = [-𝑝(1 -𝜙)]I + 2𝜇 𝑓 (1 -𝜙)D 𝑓 ( 14
)
where 𝜇 𝑓 is the asymptotic dynamic viscosity of blood. The scalar field 𝜙 is introduced to represent the volume fraction of deposited platelets (thrombus). D 𝑓 is the symmetric part of velocity gradient. The fluid phase density 𝜌 𝑓 is defined as

𝜌 𝑓 = (1 -𝜙)𝜌 𝑓 0 ( 15
)
where 𝜌 𝑓 0 is the density reference value of the fluid phase, b 𝑓 is the body force, v 𝑓 is and v 𝑇 are the velocity of the fluid and thrombus phases, respectively. 𝐶 2 is the resistance constant assuming that deposited platelets are densely compact spherical particles (2.78 𝜇m), 𝑓 (𝜙) = 𝜙(1 + 6.5𝜙) is the hindrance function. For a full description of the original thrombosis model the reader is referred to the work of Wu et al. [START_REF] Wu | Multi-Constituent Simulation of Thrombus Deposition[END_REF] and Sorensen et al. [START_REF] Sorensen | Computational simulation of platelet deposition and activation: I. Model development and properties[END_REF] .

Case of Study PFA-100®

The platelet function analyzer PFA-100® (Siemens Erlangen, Germany) is a coagulation testing device that is used to assess the primary hemostasis response [START_REF] Harrison | The PFA-100®: A potential rapid screening tool for the assessment of platelet dysfunction[END_REF] . Figure 2 shows a cross section sketch of the PFA-100® test cartridge whose main component is a bio-active membrane with a central orifice of 140 microns diameter. Whole blood is aspirated through a capillary towards the membrane, as blood flows through the central orifice of the membrane, a thrombus forms until occlusion is achieved. The membrane is coated with collagen and epinephrine or ADP to promote thrombus formation. Flow is driven by a vacuum system that applies a Δ𝑝 of 4000 Pa. Shear rates inside the orifice reach values up to 6000 s -1 promoting vWF unfolding leading to preferential thrombus formation [START_REF] Mehrabadi | A Predictive Model of High Shear Thrombus Growth[END_REF] inside the membrane orifice. When the orifice is fully occluded a Closure Time (CT) is obtained which is the quantity used to diagnose patients. The normal CTs are in the range 79-139 s for epinephrine cartridges and 61-105 s for ADP cartridges [START_REF] Harrison | The PFA-100®: A potential rapid screening tool for the assessment of platelet dysfunction[END_REF] . Prolonged CT times are observed in patients with major platelet defects (Bernard Soulier syndrome, Glanzmann's thrombasthenia, thrombocytopenia, etc.) and vWF disease, or acquired vWF disease, among others [START_REF] Harrison | The PFA-100®: A potential rapid screening tool for the assessment of platelet dysfunction[END_REF] ). In addition, CT times are affected by a variety of hematological and pharmacological factors such as hematocrit and aspirin.

Baseline Thrombosis Simulation

An axisymmetric thrombosis simulation of the PFA-100® assay is the baseline case for the uncertainty quantification analysis. Figure 3A shows the setup of the simulation. A uniform velocity profile v 𝑓 = (0, 0.127, 0) m s -1 was set at the inlet boundary. The membrane and cartridge walls were set as no-slip boundary conditions. A zero gradient velocity boundary condition was considered at the outlet. The pressure the outlet boundary was set to zero, and a zero gradient boundary condition was set for all remaining boundaries. In terms of biochemical species Table 1 lists the blood constituents concentrations that were prescribed at the inlet boundary. The platelet count values are representative of patients treated with ventricular assist devices [START_REF] Steinlechner | Platelet Dysfunction in Outpatients With Left Ventricular Assist Devices[END_REF] . At the outlet and cartridge walls a zero gradient boundary condition was applied for the species concentration. To reflect the ADP-collagen coated membrane, a reactive boundary condition was prescribed over the membrane to specify the region in which platelets are allowed to deposit. In addition, a diffusive flux of ADP = 1 × 10 7 nmol m -2 was prescribed to mimic shedding of ADP coating from the membrane. The ADP flux value was determined by best fit to clinical occlusion times reported by Steinlechner et al. [START_REF] Steinlechner | Platelet Dysfunction in Outpatients With Left Ventricular Assist Devices[END_REF] . All remaining variables were derived from the literature. The diffusion coefficient for each biochemical species were taken from [START_REF] Goodman | Computational model of device-induced thrombosis and thromboembolism[END_REF] . We employed a Newtonian constitutive model for blood with a viscosity of 3.5 cP and a density of 1050 kg m -3 . The parameter values describe along Section 2.1.1 were used in the baseline thrombosis simulation.

The computational mesh was composed of 68650 hexahedral cells with a finest resolution of Δ𝑥 = 3 𝜇m located in the membrane orifice. A dual time step strategy was used to economize computational cost. A time step of Δ𝑡 𝐶𝐷𝑅 = 1 × 10 -3 s was used to solve the species equations and thrombus growth. The flow equations are solved with a time step Δ𝑡 𝐶𝐹 𝐷 = Δ𝑡 𝐶𝐷𝑅 ∕𝑟 𝐶𝐹 𝐷 = 1×10 -8 s, ensuring the Courant-Friedrichs-Levy stability condition CFL ≤ 1, this strategy has been previously used in multi-scale platelet activation simulations [START_REF] Zhang | A predictive multiscale model for simulating flow-induced platelet activation: Correlating in silico results with in vitro results[END_REF] . A parametric study with different values of the scaling factor 𝑟 𝑓 𝑙𝑜𝑤 spanning four orders of magnitude was performed to show the robustness of the approach (See Supporting Information). Before running the thrombosis simulation, a steady state flow solution was obtained to improve stability at the first instants of the thrombosis simulation. The resulting shear rate field was computed corroborating that inside of the membrane orifice induced shear rates greater than 6000 s -1 required for vWF stretching. (See Figure 3B andC.)

Figure 4A shows the thrombus formation at 100 seconds of simulation, the blood clot is visualized through a platelet volume fraction threshold > 0.8. Thrombus is formed inside the membrane orifice within the time scales observed clinically in VAD patients (80-150 s) [START_REF] Steinlechner | Platelet Dysfunction in Outpatients With Left Ventricular Assist Devices[END_REF] , the simulation closure time is defined at 80 % of occlusion since full thrombotic closure is not numerically

Species Concentration RP 216 ×10 3 Plt 𝜇𝐿 -1 AP 2.16 ×10 3 Plt 𝜇𝐿 -1 vWF c 1000 nmol m -3 PT 1.1×10 6 nmol m -3 TB 0 nmol m -3 AT 2.844 ×10 6 nmol m -3 ADP 0 nmol m -3 TxA 2 0 nmol m -3
TABLE 1 Blood constituents inlet concentrations, the platelet count is representative of patients treated with ventricular assist devices [START_REF] Steinlechner | Platelet Dysfunction in Outpatients With Left Ventricular Assist Devices[END_REF] . achievable due to exponential increase of mechanical shear prompting thrombus cleaning and preventing full occlusion. The main mechanism that drives clot formation can be explained as follows: 1) platelets deposit on the collagen coated membrane, 2) ADP within the membrane amplifies platelet activity, 3) large shear rates inside the membrane promote both platelet activation and vWF unfolding, further amplifying platelet deposition, 4) as the clot grows, the orifice area is reduced further increasing shear rates, which increases vWF activity forming a positive feedback loop mechanism. The quantities of interest (QoIs) chosen for the UQ analysis, described below, are: (1) the percentage of occlusion over the top edge plane of the membrane and (2) the total thrombus volume in the domain. These quantities are plotted in Figure 4B.

In order to show the importance of including the vWF unfolding dynamics and its influence to the platelet adhesion process, a simulation imposing 𝑣𝑊 𝐹 𝑐 inlet condition to zero was performed. This simulation corresponds to the original model of Wu. Figur 5 A shows the absence of clot without vWF. The platelet deposition rates for both cases are presented in Fig. 5 B showcasing that vWF activity increases the deposition rates inside the membrane promoting thrombus formation.

Uncertainty Quantification

Identification and modeling of random sources

The thrombosis model used in this study is the most realistic we could assemble. In the following we will therefore put aside the model uncertainty in terms of its structure, complexity and underlying numerical algorithmic and we will focus only on physical parametric uncertainty. The model involves a total of 68 parameters (diffusion coefficients, reaction rates, biochemical concentrations, viscosity, density, etc). The uncertainties related to the microfluidic assay geometry were not considered as they are not as important as in physiological vessel hemodynamics. Due to the computational cost of performing an analysis considering all the model parameters ( 64 core-hours per simulation), we have selected a subset of 15 influential parameters based on our experience gained from previous thrombosis simulations in both macroscopic and microscopic configurations [START_REF] Rojano | Multi-constituent simulation of thrombus formation at LVAD inlet cannula connection: Importance of Virchow's triad[END_REF][START_REF] Wu | Multi-Constituent Simulation of Thrombus Deposition[END_REF][START_REF] Wu | High fidelity computational simulation of thrombus formation in Thoratec HeartMate II continuous flow ventricular assist device[END_REF][START_REF] Zhussupbekov | A Continuum Model for the Unfolding of von Willebrand Factor[END_REF] . The parameters are listed in Table 2 with their definition and distribution values. The literature does not show evidence of particular type of uncertainty distributions for these parameters. At best only low-order statistics are available (i.e. mean and standard deviation), which implies (under a maximum entropy principle) to model the distributions as normally distributed. For sake of simplicity, independent Gaussian probability distributions were assumed for the parameters. Moreover, in order to introduce a similar level of uncertainty for all input parameters, we set the standard deviation to be ten percent of the mean value for each parameter, cf. Table 2. This choice was done to avoid bias in the selection of the standard deviation value for each parameter as they are not well characterized in the literature.

Parameter

Polynomial chaos surrogate modeling framework

Numerical simulation of blood flow problems are nowadays more common, mainly because of increasing computational power, progress in imaging and numerical algorithms. Nevertheless, it remains very costly to quantify the link between multi-physics, multi-scale and multi-parameters phenomena to useful clinical indicators. A large number of model order reduction techniques have been proposed, e.g. [START_REF] Manzoni | Model reduction techniques for fast blood flow simulation in parametrized geometries[END_REF][START_REF] Adjoua | Reduced-order modeling of hemodynamics across macroscopic through mesoscopic circulation scales[END_REF] .

In order to alleviate the computational cost of a parametric study involving the full thrombosis model, a surrogate model will be constructed in place of the direct thrombosis model. The surrogate model should approximate the quantity of interest 𝑌 as accurately as possible. We remind the reader that this is not an obvious task in case one is interested in multiple quantities of interest. A continuous surrogate is constructed from a finite set of the thrombosis model predictions for some selected parameter values in a "off-line" stage, and subsequently used "in-line" during the prediction stage. To this end, different methods are available such as Gaussian processes [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] , support vector machines [START_REF] Smola | A tutorial on Support Vector Regression[END_REF] , stochastic interpolation [START_REF] Tatang | An efficient method for parametric uncertainty analysis of numerical geophysical models[END_REF] , or stochastic spectral methods such as polynomial-based representations. Polynomial chaos expansions (PCE) will be retained to express the surrogate model in a closed form [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF][START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF] and some formulations exist to adaptively construct the approximation functionals satisfying some optimality criteria, e.g. [START_REF] Lucor | Cardiovascular Modeling With Adapted Parametric Inference[END_REF] . Let (Ω, , ) be the probability space where Ω is the space of random events 𝜔, this domain has a 𝜎-algebra  and is equipped with a probability measure . The vector of normally-distributed random parameters can be written as 𝑿 ≡ 𝑿(𝜔) = (𝑋 1 , … , 𝑋 𝑑=15 ). If we consider, at a given time instant, the 𝑑-variate and second-order random variable 𝑌 ∶  𝑿 ⊆ ℝ 𝑑 → ℝ, then 𝑌 (𝑡, 𝑿) ∈ 𝓁 2 (Ω, , ), can be expressed as the following expansion [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF] :

𝑌 (𝑡, 𝑿) = ∞ ∑ 𝑗=0 𝑦 𝑗 (𝑡) ψ j (𝑿), (16)
where ψ j (𝑿) = ∏ 𝑑 𝑖=1 𝜓 (𝑖) 𝑗 (𝑋 𝑖 ) are the multivariate basis functions that form a complete basis, orthonormal with respect to the probability measure 𝜌 𝑿 of the random input, and 𝜓 (𝑖) 𝑗 are the univariate Hermite polynomial functions along the 𝑖 th dimension. In practice, PCE expansion must be truncated, and writes: 1

𝑌 (𝑡, 𝑿) ≈ ∑ 𝜸∈𝛽 𝑝 𝑦 𝜸 (𝑡) ψ 𝜸 (𝑿), (17)
where 𝑦 𝜸 (𝑡) are the evolving modal coefficients corresponding to the ψ 𝜸 basis. For instance, a choice of 𝑝 = 2 corresponds to at most quadratic Hermite polynomials in the variables 𝑿. We will restrict ourselves to tensor-product polynomial spaces ℙ 𝛽 𝑝 , where 𝛽 𝑝 is an index set of "degree" 𝑝, and where 𝑃 = dim(ℙ 𝛽 𝑝 ) ≡ #𝛽 𝑝 , will denote the cardinality of the selected polynomial space. There are different ways of constructing the approximating polynomial spaces that will impact their cardinality. More specifically, in this study we will rely on 𝑞-norm type of polynomial spaces, allowing for an hyperbolic truncation of the approximation basis (therefore reducing the computational cost):

ℙ 𝛽 𝑝,𝑞 with index set 𝛽 𝑝,𝑞 = {𝒚 ∈ ℙ 𝛽 𝑝 ∶ ‖𝒚‖ 𝑞 = ( ∑ 𝜸 𝑦 𝑞 𝜸 ) 1∕𝑞 ≤
𝑝} [START_REF] Shen | Sparse Spectral Approximations of High-Dimensional Problems Based on Hyperbolic Cross[END_REF] . By default, the hyperbolic truncation with 𝑞 = 1 reduces to the standard total-degree truncation scheme. Decreasing the value of 𝑞 de facto decreases the number of polynomials of high interaction order kept in the expansion.

Adaptively constructing the surrogate by least-squares minimization with regularization

Based on a training set {𝑿 (𝑖) , 𝑌 (𝑖) } 𝑖=1,…,𝑁 collected at each time instant, ordinary linear regression can be used to compute the unknown coefficients 𝒚 ≡ 𝑦 𝜸 , [START_REF] Choi | Polynomial Chaos expansion with Latin Hypercube Sampling for estimating response variability[END_REF][START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF] by minimizing the residuals 𝒓 ≡ 𝒀 -𝝍 𝛽 𝑝 𝒚 in the 𝓁 2 -norm through an optimization problem:

𝒚 = argmin 𝐲∈ℝ 𝑃 ‖𝒀 -𝝍 𝛽 𝑝 𝐲‖ 2 2 , ( 18
)
where 𝝍 𝛽 𝑝 is the measurement matrix corresponding to the gPC expansion in the index set 𝛽 𝑝 . The solution is obtained in matrix form as:

𝒚 = ( 𝚿 𝖳 𝛽 𝑝 𝚿 𝛽 𝑝 ) -1 𝚿 𝖳 𝛽 𝑝 𝐘, ( 19
)
where 𝐘 is a vector of observations of size 𝑁 × 1, 𝚿 𝛽 𝑝 the measurement matrix of size 𝑁 × 𝑃 with Ψ 𝑖𝑗 = ψ j (𝑿 𝑖 ), and 𝒚 the vector of coefficients of size 𝑃 × 1. Unfortunately, when the number of simulations 𝑁 < 𝑃 , the problem is not well posed under [START_REF] Nagareddy | Inflammation and thrombosis in cardiovascular disease[END_REF] Instead of indexing the expansion on a single integer amounting to the cardinality of the entire approximation space, one can also make use of a multi-index notation that is equivalent.

If 𝛽 𝑝 is an index set for multi-index 𝜸 = (𝛾 1 , … , 𝛾 𝑑 ) ∈ ℕ 𝑑 0 , then ℙ 𝛽 𝑝 ≡ span{𝜓 𝜸 | 𝛄 ∈ 𝛽 𝑝 }.

Surrogate-based global sensitivity analysis

Once the modal coefficients are computed, moments, confidence intervals, sensitivity analysis and probability density function of the solution can be readily evaluated. Global variance-based sensitivity analysis may be performed to quantify the relative importance of each (or a group of) random input parameter to the uncertainty response of the system. The Sobol' functional decomposition of 𝑌 = 𝑓 (𝑡, 𝑿) is unique and hierarchic. We have:

𝑌 = 𝑓 (𝑡, 𝑿) = ∑ 𝒔 ⊆{1,2,…𝑁 𝑑 } 𝑓 𝒔 (𝑡, 𝑿 𝒔 ), (21)
where 𝒔 is a set of integers such that 𝑿 𝒔 = (𝑋 𝑠 1 , … , 𝑋 𝑠 𝑁 ), with 𝑁 = card(𝒔) = |𝒔| and 𝑓 ∅ = 𝑓 0 = 𝔼[𝑌 ]. In this way, the variance of the response can be decomposed as [START_REF] Efron | The jackknife estimate of variance[END_REF] :

𝕍 (𝑌 ) ≡ 𝜎 2 (𝑡) = ∑ 𝒔 ⊆{1,2,…𝑁 𝑑 } 𝜎 2 𝒔 (𝑡), with 𝜎 2 𝒔 (𝑡) = 𝕍 ( 𝔼 [ 𝑌 |𝑿 𝒔 ]) - ∑ 𝒖⊂𝒔 𝒖≠𝒔 𝒖≠∅ 𝜎 2 𝑢 (𝑡), (22)
where 𝕍 and 𝔼 are the variance and the expectation operators, respectively. The normalized Sobol' indices 𝑆 𝒔 are defined as [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] :

𝑆 𝒔 ≡ 𝜎 2 𝒔 (𝑡) 𝜎 2 (𝑡)
and

∑ 𝒔 ⊆{1,2,…𝑁 𝑑 } 𝒔≠∅ 𝑆 𝒔 (𝑡) = 1, (23)
which measure the sensitivity of the variance of 𝑌 due to the interaction between the variables 𝑿 𝒔 , without taking into account the effect of the variables 𝑿 𝒖 (for 𝒖 ⊂ 𝒔 and 𝒖 ≠ 𝒔). In this work, we consider first-and total-order Sobol' indices. The firstorder quantifies the effect of the single parameter alone 𝑋 𝑠 on the output. The total index of input variable 𝑋 𝑠 , denoted 𝑆 𝑇 𝑠 , is the sum of all the Sobol' indices involving this variable: 𝑆 𝑇 𝑠 = ∑ 𝒔 ⊆{1,2,…𝑁 𝑑 } 𝑆 𝒔 . This latter definition is not practical since it would result in computing each index separately. Instead, by denoting 𝑆 s the sensitivity measure of all the variables excluding variable 𝑋 𝑠 , the total index can be rewritten as: 𝑆 𝑇 𝑠 = 1 -𝑆 s. Previous Sobol' indices can be evaluated by Monte-Carlo estimators, but require very large samples and introduce errors due to sampling. Instead, they can be obtained very straightforwardly directly from the polynomial chaos modal coefficients. Thanks to the orthonormality of the approximation basis, the Sobol' indices are explicitly obtained as combination of the PCE coefficients [START_REF] Marelli | UQLab user manual -Polynomial chaos expansions[END_REF] . Despite this "exact" relation, the Sobol' indices are only accurate if the surrogate PCE is itself accurate.

RESULTS

Multi-constituent Thrombosis Simulations

The evolution in time of closure percentage and thrombus volume is shown in Fig. 6 for 120 simulations. Closure percentage (Fig. 6 A) is observed to follow a sigmoidal curve reaching a plateau of about 90% at time of 100-200 seconds. The variance among simulations appears to be greatest during the rapid growth portion of the curve, and virtually nil during the first and last 25 seconds. It is noteworthy that the simulations do not reach full occlusion due to extremely high shear rates as the orifice cross section area approaches zero, and thrombus cleaning predominates.

In terms of the thrombus volume, Fig. 6B shows little growth during the first 100 seconds, the curves start to deviate significantly around 125 seconds reaching different levels of total thrombus volume between approximately 0.1 to 1.7 𝑚 3 × 10 -3 .

Anticipating on the results of the following settings, both plots in Fig. 6 also show the surrogate mean, std-based confidence interval, and a box plot of the 25 𝑡ℎ and 75 𝑡ℎ percentiles. It can be observed that the surrogate model is able to correctly reproduce the trend observed in multi-constituent thrombosis simulations (gray plots).

PCE surrogate model training and thorough validation

In this paper, we aimed at obtaining a subtle balance between cost and predictive accuracy. We first carried out a preliminary study to determine a very moderate number of simulations to build our training database in order to construct a metamodel with small errors on various quantities. In addition, we were very drastic in our way of validating the metamodel by several means. We relied on LOO cross-validation error. We have also pursued a generalization test on a sample that was not used during training and that includes some extreme points (data points far away from the central tendencies of our training dataset). By doing so, we have quantified the confidence in these predictions. Finally, as we will see in a later section, we have compared our predictions with clinical data, cf. Section 3.4. We will describe in the following, the details of the various analyses. To ensure that the parameter space 𝑿 was uniformly explored, a Sobol' quasi Monte Carlo sequence was used to generate a DoE matrix. This choice is based on the fact that the classic Monte Carlo method does not depend on the number of dimensions, its convergence is slow and recent papers have shown that when applied to sensitivity analyses, Sobol' sequences exhibited faster convergence. While the observed benefit to convergence is moderate, Sobol' sequences are computationally less expensive to compute than Latin Hypercube Samples (LHS) and also have the benefit of being deterministic, which allows for better reproducibility of results. In addition, a preliminary comparison was done using the aforementioned techniques (See Supplementary Information.) Then we computed the LOO error as a function of the training sample size using as QoIs the closure time and the thrombus volume as defined in Section 2.3. (See Figure 7.) This was repeated for five surrogate models of increasing complexity (increasing polynomial order 𝑝 with truncation parameter, 𝑞 = 1). Some readers may notice that some convergence curves with fixed polynomial order do not exhibit monotonous decrease despite increase in the training dataset size. This is due to the adaptive LARS approach, promoting sparsity by trimming the approximation basis. Small size basis (i.e. for low polynomial order) are sometimes affected by this optimization, resulting in small error fluctuations. Except for these very few models, the accuracy improved with increasing training size, reaching a point of diminishing returns at approximately 80-120 simulation.

In terms of model complexity, the best performance was found for polynomial order from 𝑝 = 2 and 3. Further increasing the polynomial order was found to deteriorate accuracy, reflecting the limited sample size leading to overfitting. This motivated the use of an adaptive approach that searches, for each time step, an optimal combination within given ranges, i.e. polynomial orders 𝑝 ∈ {1, 6} and truncation parameters 𝑞 ∈ {0.5 ∶ 0.1 ∶ 1}. This yielded lower error than fixed polynomial order and truncation parameter approximations. In any case, it should be reminded that with more data the surrogate model would increase in accuracy. Nevertheless, this pilot study demonstrated that 120 simulations were sufficient in order to reach reasonable accuracy when combined with an adaptive polynomial approach. When considering the closure percentage and the thrombus volume QoIs, a series of PCE meta-models were constructed, one for each time step at which the QoI was saved during the simulations. It might be useful to the reader to know that surrogates for time-dependent processes can also be constructed as in Mai and Sudret 66 , however, this is out of the scope of the current study. To cross-validate the model for each time step, 125 points were sampled using the quasi-random Sobol' sequence, of which 120 points were used for the training set to build the surrogate model and 5 points used as the test set. For these 5 points, we have made sure that the testing data were sufficiently far away from the central tendency of the joint input distribution. The Mahalanobis distance was used to quantify how far the cross-validation points are from the original training dataset (See Table 3). All the points are at least 3 stds away from the mean of the training data-set, two of them are at almost 5 stds away. Five predictions using this methodology are shown in Figure 8, the dots correspond to the thrombosis simulation, the solid line is the PCE surrogate and the shaded area is the PCE confidence interval computed using the bootstrap technique [Marelli and Sudret 67 ]. Producing a single prediction from each PCE surrogate is not very realistic due to the small training database and adaptive procedure used. Therefore, a confidence interval can be generated accordingly for each prediction. It is noteworthy that each surrogate model does not have the same selected polynomial approximation space (e.g. not same optimal order 𝑝 opt ), depending on the complex nature of the response to the parameters randomness along time. There are different approaches to combine polynomial chaos expansion and bootstrapping techniques. For the fast algorithm approach used for a large training sample, the sparse polynomial basis identified by the LARS algorithm during calibration is calculated only once from the available full experimental design, and bootstrapping is applied only to the final step, which consists in a classic ordinary least-square regression on the sparse basis. The full bootstrapping approach, recommended for smaller training sample, is more expensive but more accurate as it includes the estimation of the sparse PCE basis for each of the bootstrapped experimental designs. In this paper, the second approach is followed and prediction confidence intervals (±2std) are evaluated with 200 bootstrap replications. We notice that the simulations obtained with the full thrombus model do belong to these confidence intervals. The latter are wider for the thrombus volume variation. The accuracy of the PCE, evaluated by cross-validation with the LOO error criteria, remained for the occlusion percentage output under 20% for most cases with only one time step at 30% error (at 𝑡 = 200 seconds). In the case of the thrombus volume output the LOO error remained under 10% reaching its maximum value by the last time step of the simulation. We consider that these errors showed that our surrogate was efficient considering that our initial space has 𝑑 = 15 dimensions and we only used a moderate computational budget of 𝑁 = 120 simulations.

In conclusion, these results provided sufficient validation for our PCE approach to continue with the sensitivity analysis and further applications.

Sobol' Indices

The first order and total Sobol' indices for the 15 parameters computed using the closure percentage as QoI are shown in Fig. 9. In order to scale the effect of each parameter, the first Sobol' indices were multiplied by the variance computed through the PCE approximation. Considering the first order indices, seven parameters have a significant impact in the range 50𝑠 < 𝑡 < 200𝑠 which is in line with the fact that most variations in the QoI are observed within this time range (see Fig. 6 A). The physical significance of the parameters identified by the first Sobol' indices is related to three main factors:

• Thrombus constitution. Resting platelet count (parameter 1) is the main constituent of the thrombus, therefore, it is understandably influential in the closure percentage.

• vWF activity. The parameters 3, 8, 11, 12, and 14 (𝑣𝑊 𝐹 𝑐 , ̇𝛾𝑣𝑊 𝐹 , 𝑡 𝑣𝑊 𝐹 , 𝑊 𝑖 𝑐𝑟𝑖𝑡 𝑒𝑓 𝑓 , and 𝑣𝑊 𝐹 𝑐𝑟𝑖𝑡 𝑠 ) control the unfolding dynamics of vWF and its subsequent impact in platelet deposition. The fact that these parameters were flagged signals the important role of vWF in high shear rate thrombosis.

• Platelet-platelet adhesion. The parameter 6 (𝑘 𝑎𝑎 ) is the rate of platelet to platelet adhesion which directly influences the thrombus growth rate.

The Total Sobol' indices which evaluate the total effect of an input parameter, including all variance caused by its interactions with other parameters, are shown in Fig. 9 B. Significant effects for 11 out of 15 parameters are observed indicating nonnegligible parameter interactions. Interestingly large effects are observed specifically after 180 seconds even-though the QoI shows little variation after this time. This could be explained by the high order approximations that were retained for these time steps. Figure 10 shows the first order indices multiplied by the PCE variance and the Total Sobol' indices for the total thrombus volume QoI. In this case, only six parameters showed significant contribution to the output variations mainly at the last time steps which is in agreement with the behavior observed in the time evolution of the total volume QoI (see Fig. 6B). As in the previous QoI, the parameters that were identified are related to the thrombus constitution, the rate of platelet adhesion and vWF activity. These six parameters were also identified by the first Sobol' indices for the closure percentage QoI. Parameter 8 ( ̇𝛾𝑣𝑊 𝐹 ) that regulates vWF unfolding is not found to be influential in this case; although, a small contribution in the middle run of the thrombus formation is observed. This suggests that its role is only related to the initiation of the thrombus formation and not the subsequent growth of the thrombus. This is corroborated by the large value of its Total Sobol' index. (See Figure 10B.) The other six parameters found to be influential by the first Sobol' indices are also found to be influential considering their total Sobol' indices. Parameters 9 and 7 (Δ ̇𝛾 and 𝑘 𝑟𝑎 ) show a low but non-negligible impact, these parameters are related to direct vWF activity and platelet deposition respectively.

Clinical Application: PFA-100® Closure Time

To investigate the practical relevance and the generalization of the current UQ framework, the surrogate model predictions for three clinical scenarios were compared to reference clinical data: normal blood, thrombocytopenia, and vWF disease. The PCE surrogate model was calibrated using the closure time of the PFA-100® as QoI as defined in Section 2.3. The full set of thrombosis simulations (125) were used to compute the PCE coefficients. Three sets of 100,000 points (parameters) as defined in Table 4 were drawn using the Monte Carlo sampling technique as implemented in UQLab framework. The 13 supplementary model parameters that do not appear in Table 4 were held constant, as defined in Table 2. Once the input databases were created for each case. The surrogate model was used to predict the 100,000 closure times for each set of parameters. The probability distribution of the PCE closure times for each scenario is presented in Figure 11 (solid red lines) superimposed with histograms derived from clinical data corresponding 172 patients, reported by Harrison et al. [START_REF] Harrison | The PFA-100®: A potential rapid screening tool for the assessment of platelet dysfunction[END_REF] . It should be noted that for the vWFD and Thrombocytopenia cases, the data were consolidated, as closure times equal to 300 were not taken into account since they are the maximum limit of the coagulation test indicating that the coagulation test failed. To facilitate the comparison a distribution was computed for the clinical data: a maximum likelihood estimate of a gaussian distribution was computed for the physiological case (large enough sample) and a kernel fit was obtained for the vWFD and thrombocytopenia cases (small sample). Even though the vWFD and Thrombocytopenia cases may not be statistically significant (𝑛 = 26 and 𝑛 = 6 respectively) it can be observed that the surrogate model captures well the trend observed in the coagulation assay for the three cases; although for the thrombocytopenia case the PCE prediction lags the clinical data slightly.

Case

DISCUSSION

Substantial progress has been made over the past decade in mathematical models of thrombosis which has enabled numerical simulation of a variety of clinically relevant thrombotic scenarios [START_REF] Peach | Personalizing flow-diverter intervention for cerebral aneurysms: from computational hemodynamics to biochemical modeling[END_REF][START_REF] Ab Naim | Flow pattern analysis in type B aortic dissection patients after stent-grafting repair: Comparison between complete and incomplete false lumen thrombosis[END_REF][START_REF] Yazdani | Integrating blood cell mechanics, platelet adhesive dynamics and coagulation cascade for modelling thrombus formation in normal and diabetic blood[END_REF] . However, inherent uncertainties and assumptions of thrombosis models limit their translation to clinical practice and medical device development. This study investigated the impact of model parameter uncertainty in a thrombosis simulation by means of a polynomial chaos approximation. After extensive validation of the polynomial approximation, a global sensitivity analysis of 15 model parameters was performed. Six parameters showed significant influence in both quantities of interests considering their first Sobol' indices. This is in line with a recent study conducted by Melito et al. [START_REF] Melito | Sensitivity Analysis of a Phenomenological Thrombosis Model and Growth Rate Characterisation[END_REF] in which the authors found that only 4 of the 9 parameters in the model of Menichini [START_REF] Menichini | Predicting false lumen thrombosis in patient-specific models of aortic dissection[END_REF] were influential for predicting volume fraction of thrombus in a backwards step. However, investigation of parameter interaction via the total Sobol' indices revealed influence of a subset of the 9 "non influential" parameters according to the first Sobol' indices. This suggest that simplification of the model must consider both first and total Sobol' indices, lest it neglects the interdependence of parameters.

Several limitations are acknowledged in the current study. First, the thrombosis model itself involves several approximations [START_REF] Wu | Multi-Constituent Simulation of Thrombus Deposition[END_REF] , in which blood is treated as a Newtonian fluid, and the roles of red blood cells and fibrin formation are not explicitly considered. Yet, in reality, the closure time is likely to be influenced by RBC trapping within the clot, and possibly influenced by platelet margination induced by RBCs. Second, a first triage of parameters was perform selecting only 15 parameters to perform the sensitivity study. Even though this selection was informed by our previous thrombosis studies [START_REF] Rojano | Multi-constituent simulation of thrombus formation at LVAD inlet cannula connection: Importance of Virchow's triad[END_REF][START_REF] Wu | Multi-Constituent Simulation of Thrombus Deposition[END_REF][START_REF] Wu | High fidelity computational simulation of thrombus formation in Thoratec HeartMate II continuous flow ventricular assist device[END_REF][START_REF] Zhussupbekov | A Continuum Model for the Unfolding of von Willebrand Factor[END_REF] it remains somewhat arbitrary. Furthermore, the parameters excluded in this first triage step can have non-negligible parameter interactions with the 15 parameters selected. A possible way to better select this set of parameters is to perform a screening global sensitivity analysis with the Morris method as in Méndez Rojano et al. [START_REF] Rojano | Kinetics of the coagulation cascade including the contact activation system: sensitivity analysis and model reduction[END_REF] . Third, the type of distribution and the magnitude of its variance that was assumed for the model parameters (10%) was somewhat arbitrary and was applied to all the parameters equally. In future studies, each parameter should be assigned its own characteristic variance. In addition, the definition of closure time used in the simulations was based on 80% occlusion of the orifice, whereas in reality it would be 100%. This was necessary due to the simulations reaching an asymptotic value of approximately 90% due to extremely high shear rates as the orifice diameter approaches zero. Finally, the clinical database for normal blood was constructed with n=172 patients; for the vWFD and thrombocytopenia, only n=26 patients and n=6 patients were available respectively. A larger sample size would improve the comparisons for vWFD and thrombocytopenia cases.

Despite its limitations, the PCE surrogate model demonstrated practical utility inasmuch as it permitted quantification of the uncertainty associated with a large number of parameters of a thrombosis model with economical computational cost. By contrast, it would have been very prohibitive to simulate all 300,000 sets of parameters with the full thrombosis model. Moreover, it was shown that the adaptive surrogate initially calibrated for certain parametric variabilities was nevertheless also capable of providing very reasonable predictions in parametric regions of lower probabilities. It is particularly noteworthy that despite the fact that the parametric distribution of the normal clinical case was centered outside the training confidence interval (estimated for this normal distribution as 1.3 times the RP mean training value = 280𝑒3 Plt 𝜇 L -1 ), the surrogate model nevertheless predicts very well the clinical observations.

In conclusion this study showed that a PCE surrogate model trained with high fidelity thrombosis simulation can reproduce the trends observed in clinical data. These results encourage future studies using this methodology to investigate the influence of anti-platelet agents to identify scenarios that could achieve a balance between the risk of thrombosis and bleeding.

ACKNOWLEDGMENTS

This research was supported by NIH R01 HL089456 and NIH R01 HL086918.

Bibliography

FIGURE 1 A

 1 FIGURE 1 A) Diagram of platelet activation and deposition in the thrombosis model. 𝑘 𝑟𝑝𝑑,𝑏 and 𝑘 𝑎𝑝𝑑,𝑏 are the rates at which resting and activated platelets deposit to the surface. Resting platelets can be activated by mechanical shear, or the combination of agonists: ADP, TxA 2 and thrombin. Once deposited, AP 𝑑 can stabilize or embolize due to flow shearing forces. B) Schematic depiction of vWF unfolding and vWF-mediated platelet deposition and aggregation. The presence of stretched vWF s amplifies the deposition rate of free-flowing platelets by 𝛼 𝑘 and increases the resistance of deposited platelets to shear cleaning by 𝛼 𝜏 .

FIGURE 2

 2 FIGURE 2 Cross section of PFA-100® testing cartridge showing the capillary, the central membrane with orifice, and reservoir. As whole blood is aspired through the cartridge blood constituents aggregate in the coated membrane orifice. The dimensions of the PFA-100 cartridge as considered in the simulations are: 𝐷 𝑐 = 200, 𝐷 𝑚 = 147,𝐷 𝑟 = 1200, 𝐻 1 = 1000, 𝐻 2 = 400, 𝐻 3 = 1200, and 𝐻 𝑚 = 200 all in 𝜇m.

FIGURE 3 A

 3 FIGURE 3 A) Numerical setup of PFA-100® thrombosis simulation. The membrane is set as reactive boundary condition allowing platelet deposition. B) Scalar shear rate field, large shear rates are observed in the membrane orifice. C) concentration of stretched vWF (vWFs) illustrating unfolding within and downstream of the membrane orifice.

FIGURE 4 A

 4 FIGURE 4 A) 3D rotational extrusion view of thrombus formation at 100 seconds of simulation, clot forms inside the membrane orifice. The blue line and arrow denote the rotation axis to obtain the 3D view. B) Quantities of interest (𝑌 (𝑋 𝑖 )) for the UQ analysis, occlusion percentage across the top edge of the membrane (dashed line) and total thrombus volume inside (solid line).

FIGURE 5

 5 FIGURE 5 Thrombosis simulation with inlet condition 𝑣𝑊 𝐹 𝑐 = 0 which corresponds to the original model of Wu et al. 33 . A) Lack of thrombus formation at 100 seconds whith inlet boundary condition 𝑣𝑊 𝐹 𝑐 = 0. B) Platelet deposition rate 𝑘 𝑎𝑝𝑑 with and without 𝑣𝑊 𝐹 𝑐 .

FIGURE 6

 6 FIGURE 6 Output variables for 120 thrombosis simulations (gray lines), the blue solid and dashed lines are the mean and standard deviation of the surrogate model respectively, the box plots show the median and 25 𝑡ℎ and 75 𝑡ℎ percentiles of the surrogate model predictions. A) Occlusion percentage evolution in time B) Thrombus volume in the computational domain in time.

FIGURE 7

 7 FIGURE 7 Time-averaged leave-one-out (LOO) cross-validation error of the parametric models as a function of the training sample size and model complexity. A) closure percentage error, and B) Thrombus Volume error. Adaptive metamodeling construction are more efficient than models with fixed complexity. They show reasonable accuracy for a budget of 120 full thrombosis simulations.

FIGURE 8

 8 FIGURE 8 Cross validation on extreme events: generalization from PCE surrogate models (mean: triangles) compared to thrombosis simulations not included in the training dataset (circles). The shaded areas, obtained thanks to bootstrapping technique, represent the (2-std) confidence intervals associated with the meta-models. Most of the thrombosis simulations fall within the confidence intervals.

FIGURE 9

 9 FIGURE 9 Sobol' indices evolution over time for the QoI closure percentage. A) First order indices multiplied by the output PCE variance. B) Total Sobol' indices.

FIGURE 10

 10 FIGURE 10 Sobol' indices evolution over time for the thrombus volume QoI. A) First order indices. B) shows the total Sobol' indices.

  count  (300𝑒3, 30𝑒3) Plt 𝜇𝐿 -1 vWFD vWF c (70 % of vWF concentration)  (700, 70) nmol m -3 Thrombocytopenia RP (50% of Normal Resting platelet count)  (150𝑒3, 15.0𝑒3) Plt 𝜇𝐿 -1

FIGURE 11

 11 FIGURE11 Distribution of closure times for simulated and clinically reported PFA-100® A) physiological case, B) vWF disease, and C) thrombocytopenia. Black circles along horizontal axes are the clinical data extracted from[START_REF] Harrison | The PFA-100®: A potential rapid screening tool for the assessment of platelet dysfunction[END_REF] , a blue histogram was constructed using this data. The red line is the probability density function approximated from one hundred thousand realizations of the PCE surrogate model, the blue dashed line in A shows a Gaussian approximation (with a maximum likelihood estimation), and the green dashed lines in B and C are the kernel fits.

  If 𝑊 𝑖 𝑒𝑓𝑓 falls below 𝑊 𝑖 𝑒𝑓𝑓,𝑐𝑟𝑖𝑡 , the stretched vWF s exhibits hysteresis and does not collapse back to vWF c until 𝑊 𝑖 𝑒𝑓𝑓 < 𝑊 𝑖 𝑒𝑓𝑓,ℎ𝑦𝑠𝑡 , as shown in Eq.(11).

						, while
	𝑘 𝑠𝑡𝑒𝑎𝑑𝑦 𝑠-𝑐	is zero.
				𝑊 𝑖 𝑒𝑓𝑓 =	√ 𝜆 (‖𝑺‖ + ‖𝜴‖) 𝜏 𝑟𝑒𝑙	(9)
		𝑘 𝑠𝑡𝑒𝑎𝑑𝑦 𝑐-𝑠	=	{ 𝑘 ′ 𝑊 𝑖 𝑒𝑓𝑓 𝑊 𝑖 𝑒𝑓𝑓,𝑐𝑟𝑖𝑡 𝑘 𝑠ℎ𝑒𝑎𝑟 𝑐-𝑠 ,	, if 𝑊 𝑖 𝑒𝑓𝑓 ≥ 𝑊 𝑖 𝑒𝑓𝑓,𝑐𝑟𝑖𝑡 if 𝑊 𝑖 𝑒𝑓𝑓 < 𝑊 𝑖 𝑒𝑓𝑓,𝑐𝑟𝑖𝑡	(10)
		𝑘 𝑠𝑡𝑒𝑎𝑑𝑦 𝑠-𝑐	=	{ 0, if 𝑊 𝑖 𝑒𝑓𝑓,ℎ𝑦𝑠𝑡 ≤ 𝑊 𝑖 𝑒𝑓𝑓 < 𝑊 𝑖 𝑒𝑓𝑓,𝑐𝑟𝑖𝑡 . 𝑘 ′ , if 𝑊 𝑖 𝑒𝑓𝑓 < 𝑊 𝑖 𝑒𝑓𝑓,ℎ𝑦𝑠𝑡 .	(

TABLE 2

 2

	𝑋 𝑖	Definition	Distribution (mean, variance)
	1) RP	Resting platelets	 (216𝑒3, 21.6𝑒3) Plt 𝜇𝐿 -1
	2) AP	Activated platelets	 (2.16𝑒3, 0.216𝑒3) Plt 𝜇𝐿 -1
	3) vWF c	vWF concentration	 (1𝑒3, 0.1𝑒3) nmol m -3
	4) 𝐾 𝑟𝑝𝑑	Deposition rate RP (to membrane)	 (1.0𝑒-8, 1.0𝑒-9) m s -1
	5) 𝐾 𝑎𝑝𝑑	Deposition rate AP (to membrane)	 (3.0𝑒-6, 3𝑒-7) m s -1
	6) 𝐾 𝑎𝑎	Deposition rate of AP to 𝐴𝑃 𝑑	 (3.0𝑒-5, 3𝑒-6) m s -1
	7) 𝐾 𝑟𝑎	Deposition rate of RP to 𝐴𝑃 𝑑	 (3.0𝑒-6, 3𝑒-7) m s -1
	8) ̇𝛾𝑣𝑊𝐹	critical shear rate for vWF unfolding	 (5.5𝑒3, 0.55𝑒3) s -1
	9) Δ ̇𝛾	vWF shear unfolding transition width	 (1.271𝑒3, 0.1271𝑒3) s 2
	10) 𝜆 𝑐𝑟𝑖𝑡	Flow type critical value for strong unfolding	 (5𝑒-3, 5𝑒-4) dimensionless
	11) 𝑡 𝑣𝑊 𝐹	vWF relaxation time	 (0.05, 0.005) s
	12) 𝑊 𝑖 𝑒𝑓𝑓,𝑐𝑟𝑖𝑡	𝑊 𝑖 𝑒𝑓𝑓 critical value for strong unfolding	 (0.316, 0.0316) dimensionless
	13) 𝑊 𝑖 𝑒𝑓𝑓,ℎ𝑦𝑠𝑡	𝑊 𝑖 𝑒𝑓𝑓 hysteresis value	 (0.053, 0.0053) dimensionless
	14) vWF crit s 15) 𝐷 ℎ 𝑒 ∕𝐷 𝑙 𝑒	Critical vWF s concentration value Ratio of adhesive forces at high and low shear rates	 (50, 5)  (500, 50) dimensionless

Fifteen (15)

parameters of the thrombosis model.

TABLE 3

 3 Mahalanobis Distance (MD) between the cross validation points and the training dataset mean in terms of standard deviations.

	Point 121 122 123 124 125
	MD	5.32 3.97 4.87 3.3 3.86

TABLE 4

 4 Definition of random parameters for the three cases.

this formulation and one needs to add a regularization term in the form of:

with

, which is moreover a good strategy to favour sparsity of the surrogate in high dimension, i.e. that it forces the minimization to favour low-rank solutions. In this paper, we rely on one of the several algorithms used to solve this constrained optimization problem, namely least angle regression (LAR) [START_REF] Efron | Least angle regression[END_REF] , which has been adapted to the context of PCE by Blatman and Sudret [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] . They have designed an adaptive hybrid form of LAR in order to obtain a sparse PCE representation that remains reasonably accurate even with for small experimental designs. Other more robust sparse representation approaches favor the early detection of data outliers [START_REF] Van Langenhove | Robust uncertainty quantification using preconditioned least-squares polynomial approximations with 𝐿 1 -regularization[END_REF] . In practice the adaptive PCE-LAR algorithm will be run at each time instant at which the data sample from the full thrombosis model is available. For each surrogate build, it will test through a polynomial order range and will select the optimal polynomial based on some validation criteria described hereafter. Thanks to the formulation, the retained polynomial basis will most likely be quite sparse (which means that a large portion of the modal coefficients will be null), resulting in a compact surrogate model. The leave-one-out (LOO) cross-validation error is a statistical learning technique, that is conveniently used for cases with small design of experiments (DoE). It consists in building 𝑁 metamodels, each one created on a reduced experimental design 𝑿 𝑙 , i.e. the training set obtained by removing a point (𝑿 (𝑙) , 𝑌 (𝑙) ) 𝑙∈{1,2,…,𝑁} , and comparing its prediction on the excluded scenario 𝑿 (𝑙) with the real value obtained from the thrombosis simulation. These 𝑁 errors are then averaged into a single value which may be collected for several time instants and different approximation bases for instance.


Slim Zidi
email: slim.zidi02@etud.univ-paris8.fr

Nadia Hamani

Lyes Kermad

Reconfigurable Supply Chain Performance: A Bibliometric Analysis

Keywords: Reconfigurable Supply Chain, Performance Evaluation, Supply Chain Disruption, Bibliometric Analysis

Market disruptions and changes make the assessment process of supply chain performance more complex. Hence, the implementation of key indicators to evaluate the ability of the supply chain to cope with these disruptions has received increasing attention. Reconfigurability provide the ability to quickly change supply chain structure and functions by increasing its responsiveness and flexibility to deal with disruptions. This paper aims to analyze the contributions of researchers in the field of performance evaluation under supply chain disruption. A bibliometric analysis is conducted to identify publications dealing with performance evaluation under supply chain disruption field. The analysis started with the selection of articles using the following keywords (Reconfigurability, Evaluation, Supply Chain). Then, the selected articles were analyzed based on the criteria of the bibliometric analysis. The paper aims to capture the current research and potential research related to the evaluation of performance in reconfigurable supply chain under disruption.

Introduction

The need to adapt to market changes requires the implementation of new supply chain redesign strategies. Market disruptions and changing customer requirements are causes that drive managers and decision makers to restructure their supply chains. Reconfiguration of supply chains presents an effective solution to meet new market needs [START_REF] Zidi | Proposition d'une approche de reconfiguration des processus métiers de la chaîne logistique[END_REF][START_REF] Zidi | A reconfiguration approach for a supply chain tracking platform[END_REF]. The reconfigurable supply chain is a flexible chain able to change its structure with minimum resources [START_REF] Chandra | Supply Chain Configuration: Concepts, Solutions, and Applications[END_REF]. It refers to changing the structure and functionality of the supply chain in a cost-effective, responsive, sustainable and resilient manner [START_REF] Dolgui | Reconfigurable supply chain: the Xnetwork[END_REF]. The reconfigurability assessment is an important phase to measure the ability of the supply chain to cope with disruptions and meet new market requirements [START_REF] Zidi | New metrics for measuring supply chain reconfigurability[END_REF].

Several concepts have been discussed in the literature that refer to the ability of the supply chain to adjust its structure and functions with new market needs such as agility, flexibility, etc. [START_REF] Chandra | Supply Chain Configuration: Concepts, Solutions, and Applications[END_REF] consider agility and flexibility as advantages of implementing a reconfigurable supply chain. Reconfigurability is characterized by modularity, convertibility, integrability, diagnosability, scalability and customization that reduce the effort of reconfiguration. These characteristics can be considered as performance indicators to evaluate the reconfigurability of supply chains [START_REF] Zidi | New metrics for measuring supply chain reconfigurability[END_REF]. Nevertheless, the lack of works elaborated in the area of reconfigurability evaluation in the network/supply chain level leads us to study the different aspects and concepts that can be addressed in the process of reconfigurable supply chain evaluation. The aim of this paper is to study the concepts addressed in the context of the evaluation of reconfigurable supply chains and to analyze the differences between them and to respond to the following question:

 What are the criteria that can be considered to evaluate supply chain performance in the context of reconfiguration?  How can the supply chain's ability to cope with disruptions be assessed? Answering this question provides a way for managers and decision makers to have visibility into the ability of their supply chain to resist and cope with disruptions. Indeed, it requires an in-depth study of the performance indicators that can be used to better evaluate the performance of the supply chain in the context of reconfiguration. For this purpose, a bibliometric analysis has been conducted in order to gather the concepts used in the context of the above problematic and to identify the most relevant performance indicators that provide managers the possibility to evaluate and improve the performance of their supply chain in a context of disruption.

The rest of the paper is organized as follows. In Section 2, we present a definition of the reconfigurability evaluation. The research methodology is described in Section 3. Section 4 presents the analysis and the obtained results of the bibliometric analysis. Section 5 concludes the paper.

Reconfigurability Evaluation

The evaluation of the ability of the supply chain to change its current configuration by reducing the effort of reconfiguration is called the "Reconfigurability Evaluation". Reconfiguration can be a strategy for improving supply chain performance [START_REF] Ma | Modeling supply chain's reconfigurability using fuzzy logic[END_REF]. The concept of supply chain reconfiguration is defined as the structural and functional change of the supply chain. It concerns the change of supply chain configuration that is considered as a set of nodes and connections. The need for reconfiguration can be triggered by a disruptive event that forces decision makers to quickly react to ensure the functioning of the supply chain, or by a decision made by decision makers to ensure the continuous improvement of their supply chain. The need for reconfiguration is generally due to a hazard or disruptive event that causes a supply chain failure. Then the reconfiguration can be deployed within the innovation strategy to improve the operational and organizational performance of the supply chain without being affected by an external event. Reconfiguration is linked to several levels of application (network, system, plant and machine) [START_REF] Napoleone | Core Characteristics of Reconfigurability and their Influencing Elements[END_REF][START_REF] Wiendahl | Justifying Changeability. A Methodical Approach to Achieving Cost Effectiveness[END_REF]. Indeed, each level of application has its own parameters and attributes allowing to judge their degree of reconfigurability and to choose the best configuration that changes its structure and functions easily and quickly. The reconfigurability characteristics (modularity, integrability, convertibility, diagnosability, scalability and customization) can be Reconfigurable Supply Chain Performance 159 relevant indicators for the evaluation of supply chain reconfigurability [START_REF] Zidi | New metrics for measuring supply chain reconfigurability[END_REF], thanks to their effective roles in reducing the reconfiguration effort [START_REF] Koren | Design of reconfigurable manufacturing systems[END_REF][START_REF] Napoleone | A framework to manage reconfigurability in manufacturing[END_REF]. Several indicators have been used to measure the degree of reconfigurability in the different previously mentioned levels such as lead time, reconfiguration time and cost, reliability, productivity, etc. Several concepts can describe this ability, such as agility, flexibility, adaptability, alignment, etc. The performance of the supply chain is significantly dependent on the ability to adapt to the dynamic environment [START_REF] Vanderhaeghe | How to Fail at Flexibility[END_REF]. These can be indicators to measure the ability of the supply chain to cope with market disruptions.

Research Methodology

Bibliometric analysis is defined as a research collection technique that studies a specific research area quantitatively using mathematical and statistical methods [START_REF] Bellis | Bibliometrics and Citation Analysis: From the Science Citation Index to Cybermetrics[END_REF][START_REF] Ding | Bibliometric cartography of information retrieval research by using co-word analysis[END_REF][START_REF] Godin | On the origins of bibliometrics[END_REF][START_REF] Merigó | Research in Production and Operations Management: A University-Based Bibliometric Analysis[END_REF][START_REF] Pourkhani | The impact of social media in business growth and performance: A scientometrics analysis[END_REF][START_REF] Pritchard | Statistical Bibliography or Bibliometrics?[END_REF]. Also, it is the quantification of bibliographic information for use in analysis [START_REF] Garfield | Citation Data as Science Indicators[END_REF]. It is used particularly in analyzing publications' content and network [START_REF] Ellegaard | The bibliometric analysis of scholarly production: How great is the impact?[END_REF]. The bibliometric has advantages in predicting future trends of disciplines [START_REF] Wang | Integrating blockchain technology into the energy sectorfrom theory of blockchain to research and application of energy blockchain[END_REF]. The Scopus database is used for our literature search. As indicated in [21], it is "the largest abstract and citation database of peer-reviewed literature: scientific journals, books and conference proceedings". First, a set of keywords was used to determine the application domains of the reconfigurability concept and to select the words indicating the reconfigurability evaluation indicators. Secondly, this search was restricted by using the "title, abstract, keywords" search in Scopus database to determine the most important aspects related to the evaluation of the supply chain reconfigurability. The initial search yielded 2325 articles. The minimum number of occurrences of a keyword was set at 10. Only 88 of the 6500 keywords reached this threshold. The keywords "Reconfigurable", "Supply Chain", "Evaluation" and "Performance" were used to obtain a first list of articles using them. The second search was more precise and found 1149 papers using "Supply chain", "Reconfigurability" and "Criteria" as keywords. Only 101 keywords reached the 5 keyword threshold out of a total of 3324. The obtained results are summarized in Table 1. The VOSviewer software is used to visualize and explore maps based on network data obtained. The steps of the proposed approach are presented in Fig. 1. The first research aims to identify the concepts used in the literature in the context of supply chain performance evaluation. This analysis allowed to identify 6 clusters that indicate the 5 aspects referring to the ability of the supply chain to resist and cope with disruptions. Then, by including the keyword " criteria ", our first search was consistent with the results of the second search, which allowed to identify 5 relevant indicators to evaluate the performance of the supply chain in the context of reconfiguration. The first search results generated six clusters as shown in Fig. 2. This clustering groups all the selected keywords into clusters that designate mainly the domain of application and the aspects related to reconfigurability. The clusters in red and light blue concern the first aspect called "changeability" which is related to product development and sustainability. The second aspect called " Resilience & Robustness " relates to the cluster in green. The cluster in blue represents the third aspect called "Reliability". The fourth aspect is deduced from the yellow cluster and is called "Flexibility". The last cluster in purple is the fifth aspect called "Agility".

Fig. 2. Co-occurrences of keywords of evaluated papers from the first search

The changeability

The Changeability aspect is related to the product development at the level of the manufacturing system. Indeed, production is an important enabler to ensure a changeable supply chain. The implementation of a changeable supply chain requires the coordination of production capacities [START_REF] Niemann | Development of a Method to Increase Flexibility and Changeability of Supply Contracts in the Automotive Industry[END_REF]. The dynamism of the production leads to the necessity of a changeable organization and the differentiation of the products to Reconfigurable Supply Chain Performance 161 ensure the sustainability [START_REF] Wiendahl | Changeable Manufacturing -Classification, Design and Operation[END_REF]. The implementation of a changeable system provides flexibility and reconfigurability [START_REF] Francalanza | Deriving a Systematic Approach to Changeable Manufacturing System Design[END_REF]. In fact, the changeability allows to evaluate the ability to change the production capacity to quickly and economically adjust to the new requirements.

The Resilience & Robustness

The Resilience & Robustness aspect designates the ability of the supply chain to resist and avoid change and succeed after failure. [START_REF] Durach | Antecedents and Dimensions of Supply Chain Robustness: A Systematic Literature Review[END_REF] consider that resilience and avoidance are the two dimensions that allow to judge the degree of robustness of the supply chain. The ability of the supply chain to succeed after a failure refers to resilience [START_REF] Zhang | Toward a resilient manufacturing system[END_REF], which means that the supply chain is able to change and reconfigure its structure and resources to achieve the expected functions [START_REF] Zhang | On the principle of design of resilient systemsapplication to enterprise information systems[END_REF]. Resilience & Robustness can be created to mitigate threats to organizational performance caused by market disruptions [START_REF] Brandon-Jones | A Contingent Resource-Based Perspective of Supply Chain Resilience and Robustness[END_REF]. From a reconfiguration perspective, Resilience and Robustness are important indicators because they allow judging the ability of the supply chain to withstand or avoid failures (before the failure occurs) or to return to its normal state with minimal delay (after a failure occurs).

The reliability

The Reliability aspect refers to the ability of the supply chain to meet the needs of customers within the required conditions. Reliability is a key indicator for quantifying risks and uncertainties in the supply chain [START_REF] Ha | A mathematical definition and basic structures for supply chain reliability: A procurement capability perspective[END_REF]. It unifies the fuzzy and random meaning in a dynamic environment [START_REF] Miao | The uncertainty evaluation method of supply chain reliability[END_REF]. [START_REF] Mittal | An Overview of Performance Measures in Reconfigurable Manufacturing System[END_REF] consider reliability as a key indicator in the performance evaluation of reconfigurable systems as it has a positive impact on responsiveness. In the process of evaluating the reconfigurability of supply chains, reliability increases the visibility of the supply chain as well as the quality of service with customers. The implementation of new technologies, especially those of digitalization, allows to improve the reliability of the supply chain, especially the reliability related to the information flows.

The Flexibility

The flexibility aspect designates the ability of the supply chain to quickly respond to changing customer requirements. It is considered as an advantage of reconfigurable supply chain [START_REF] Chandra | Supply Chain Configuration: Concepts, Solutions, and Applications[END_REF]. Flexibility allows to deal with market uncertainty [START_REF] Zhong | Flexibility assessment and management in supply chain : a new framework and applications[END_REF]. The flexible supply chain is defined as the ability to correctly and quickly react to changes by restructuring the chain cost-effectively to cope with market disruptions and uncertainties [START_REF] Kumar | Implementation and management framework for supply chain flexibility[END_REF][START_REF] Wadhwa | Supply Chain as a Flexible System: A KM Focused Competence[END_REF][START_REF] Tachizawa | Supply flexibility strategies in Spanish firms: Results from a survey[END_REF]. In reconfigurable supply chains, as flexibility can be considered as an advantage of reconfigurability, it is necessary to measure the ability of the supply chain to respond to new customer requirements in an uncertain environment to improve reconfigurability.

The agility

The agility aspect is closely related to the flexibility aspect. This aspect is defined as the ability to quickly react to short-term changes in demand or supply [START_REF] Lee | The Triple-A Supply Chain, the-triple-a-supply-chain[END_REF].

Responding quickly to new requirements requires the integration and sharing of information within collaboration strategies between all actors in the supply chain [37][START_REF] Whitten | Triple-A supply chain performance[END_REF][START_REF] Jr | The impact of timely information on organisational performance in a supply chain[END_REF]. Agility is seen as a benefit of the reconfigurable supply chain through new technologies that allow engagement with innovative supply chain partners [START_REF] Chandra | Supply Chain Configuration: Concepts, Solutions, and Applications[END_REF]. [37] consider that without flexibility, the supply chain cannot be agile. As part of reconfiguration, the measurement of agility allows to assess the ability of the supply chain to survive with new challenges and engage in new market opportunities.

Discussion

Evaluating the performance of a reconfigurable supply chain essentially requires measuring its degree of reconfigurability. This measurement consists of assessing the ability of the supply chain to change its structure and functions to cope with market disruptions and meet new customer requirements. The results of the bibliometric analysis showed several elaborated concepts that designate the supply chain's ability to adapt to new market needs (Changeability, Resilience & Robustness, Reliability, Flexibility and Agility). Furthermore, a more precise research was developed to verify the concepts related to the evaluation of the reconfigurability of supply chains as explained in the previous section as shown in Fig 3. The results of this research gave the same concepts found in the first research. This shows the importance of the aspects previously mentioned in the evaluation process Reconfigurable Supply Chain Performance 163 of the reconfigurable supply chain by integrating the different levels such as the manufacturing system, the machine, the product, etc. In the last few years, these concepts have been addressed in the context of new technologies such as industry 4.0, digitalization, cyber physical systems, machine learning, artificial intelligence. The 5 aspects found also showed the importance of sustainability in improving the reconfigurability of supply chains.

The assessment of the supply chain's ability to cope with disruptions and change its configuration to adapt to new changes requires the consideration of indicators to judge this ability. The main changes are related to the implementation of new technologies. Although our analysis shows the relevance of the 5 aspects found, it will be interesting to study how these aspects can be ensured and improved by integrating Industry 4.0 technologies. Indeed, the parameters for evaluating the supply chain's ability to change must include the characteristics of Industry 4.0 and digitalization, namely the reduction of truck time, the time between the need for reconfiguration and its realization, the increase in visibility of information, promoting adaptation with new changes, etc. On the other hand, it is recommended to integrate the concept of sustainability and its aspects in the parameters of performance evaluation of the supply chain under disruption. Indeed, the improvement of changeability, resilience & robustness, reliability, flexibility and agility must include other aspects, in this case the aspects of Industry 4.0 and its positive impact on sustainability while promoting the possibility of easy reconfiguration.

Conclusion

In this article, a study of the aspects and concepts treated within the context of the evaluation of reconfigurable supply chains was elaborated. Based on a bibliometric study, an analysis has been carried out through two searches of papers dealing with the evaluation of reconfigurable supply chains.

The first study gave the results of five essential aspects related to the evaluation of reconfigurability (Changeability, Resilience & Robustness, Reliability, Flexibility and Agility). These aspects allow an effective assessment of the supply chain's ability to cope with disruptions and meet new market requirements. The second, more precise research validated the five aspects identified above. This study showed the crucial role and the high dependence between these aspects and Industry 4.0 technologies in the supply chain reconfiguration strategy, as it is the most common keyword found in the researches elaborated. Thus, the latter must integrate the concept of sustainability to meet its requirements. Thus, the guarantee of the six characteristics of reconfigurability (modularity, convertibility, integrability, diagnosability, scalability and customization) provides a reconfigurable supply chain and has a crucial role in ensuring changeability, resilience & robustness, reliability, flexibility and agility.

The proposed approach has two main limitations. First, the lack of detail on the parameters and metrics allowing to quantitatively measure the five aspects previously mentioned. Secondly, it is recommended to consider the interrelationships between

Fig. 1 .

 1 Fig. 1. Research methodology

Fig. 3 .

 3 Fig. 3. Co-occurrences of keywords of evaluated papers from the restricted search

these five aspects in order to better optimize the model of evaluation and improvement of the supply chain's ability to cope with disruptions.

As a perspective, we propose to study the evaluation of reconfigurable supply chains in the area of Industry 4.0.