An adjoint‐based method for the computation of gradients in coagulation schemes

An adjoint‐based methodology is proposed to compute the gradient of the outcomes of mathematical models for the coagulation cascade. The method is first exposed and validated by considering a simple, analytically tractable case involving only 3 species. Its potential is further illustrated by considering a detailed model for the extrinsic pathway involving 34 chemical species interacting through 45 chemical reactions and for which the gradient of Endogeneous Thrombin Potential, clotting time, maximum rate and peak value of thrombin with respect to the initial concentrations and reactions rates are computed. It is shown that the method produces gradients estimates that are fully consistent with the finite differences approximation used so far in the literature, but at a much lower computational cost.


| INTRODUCTION
Haemostasis is the set of phenomena at play in blood and its surrounding tissues to prevent or allow the cessation of blood flow in the vessels. 1 Disorders of this mechanism may lead to either excess bleeding or thrombus formation, thus threatening the life of the patient due to haemorrhage or thromboembolic events, respectively. Platelets play a central role in the clotting process, being activated and recruited to form the thrombus which is consolidated thanks to fibrin strands. Fibrin is produced from the action of thrombin onto fibrinogen, a protein diluted in plasma. Thrombin, an enzym (serine protease), is the end product of the coagulation cascade, a series of biochemical reactions with autocatalytic amplification and several inhibition routes. 2 Thrombin has a central role in haemostasis since, on top of cleaving fibrinogen into fibrin, it enters in several positive (pro-coagulant) feedbacks (activation of coagulation factors XIII, XI, VIII, V and platelets). Once bound to thrombomodulin, thrombin also has a negative (anti-coagulant) action through the activation of protein C which can deactivate the activated factors VIII and V. Because it has such a central role in haemostasis, thrombin is involved in different tests performed in the clinical routine to assess the thrombogenic profile of patients, such as the thrombin time assay or the thrombin generation time assay. In the latter, the time evolution of thrombin concentration in a sample of (platelet rich) plasma is measured after the adjunction of tissue factor and calcium to trigger the coagulation cascade. 3 Because it involves many biochemical species interacting through many positive/negative feedback loops, the coagulation cascade is by itself a very complex phenomenon. Mathematical models describing the kinetics between the different species are now available, as the outcome of 40 years of intense research. If these models will most probably never replace biological tests, 4 they can help improving our understanding of the coagulation cascade. 5 A comprehensive model of the so-called extrinsic pathway of coagulation has been proposed by Hockin et al. 6 and was further completed by Butenas et al. 7 The activation of the coagulation by tissue factor, like in the extrinsic pathway, is not the only route to the generation of thrombin; the intrinsic pathway, initiated by the activation of the Hageman factor (XII) on a negative external surface was modelled by Chatterjee et al. 8 In this work, the fact that many of the chemical reactions are heterogeneous (they take place efficiently only once the reactants are bound to a phospholipid surface such as the membrane of an activated platelet 9 ) is accounted for heuristically by introducing a global activation state variable upon which the values of the constant rates depend. A more physically sound approach was introduced by Fogelson et al., 10 where the number of activated platelets and available bound sites are part of the unknowns being tracked by the model, on top of the bulk and surface concentrations of the reactive species. The transport phenomena (convection of platelets, diffusion of chemicals) are not neglected in 10 but they are represented by simple laws based on global transport coefficients. 11,12 In doing so, the coagulation cascade is represented by a set of coupled non-linear differential equations stemming, in most of the cases, from the law of mass action. The resulting model reproduces many interactions relevant to the coagulation cascade and was for example able to explain how, in the case where factor V is low, the generation of thrombin can reach normal values despite a strong deficiency in factor VIII (haemophilia A). 13 The mechanism at play results from the competitive access to the binding sites on the membrane of activated platelets and was confirmed experimentally after being discovered numerically. 13 Except when the transport phenomena are explicitly accounted for and the concentrations of the species depend both in time and space (e.g., Ref. [14][15][16], mathematical models for coagulation consist in a set of ordinary differential equations whose parameters (constant rates) and initial conditions (species concentrations) are uncertain due to physiological variability, limited knowledge of some mechanisms and measurements inaccuracy. According to Danforth et al., 17 the normal ranges of coagulation factors were determined in 2006/2007 by The Fletcher Allen Health Care Special Coagulation/Haematology Laboratory (Burlington, Vermont) from 150 volunteers; if assessed as the mean AE twice the standard deviation, these ranges are of order a few tens of percent, AE40%, say. Such a assessment does not seem to be available for the constant rates and variations in the range 10% À 1000% were considered in Danforth et al., 18 whilst 50% À 150% were used in Link et al., 19 without further justification.
Making the models useful in the current efforts towards a better understanding of coagulation requires knowing how thrombin production depends on the uncertain inputs. The Monte Carlo approach, where about 10 4 À 10 5 simulations are performed with randomly selected input values, was followed by a few groups in the past. 17,18,[20][21][22] This methodology is accurate but time consuming and hardly tractable for complex models. Several parametric studies were also performed, where the outcome of the model is analysed for different deficiencies in coagulation factors. 6,8,10,23,24 Again, this requires computing many different situations, the number of which increases dramatically with the number of parameters in the model.
The gradient of the outcome contains relevant information about the change in the model output with respect to any change in the input parameters (either constant rates or initial concentrations). This vector contains the basic information required to perform a local, one-at-a-time (OAT) sensitivity analysis. Contrary to a global sensitivity analysis which considers a whole range of operating points and simultaneous variations of input variables, a local approach examines how the system evolves around a baseline case when its parameters are modified by disturbances of low amplitude; the OAT approach analyses how the outcome changes when only one parameter departs from its baseline value. Although fully justified for linear problems only, 25 a local OAT analysis is often performed, either as a first step before a more complete and CPU demanding global analysis, or as a final product to obtain useful, albeit incomplete, information on sensitivity of the system of interest. For example, in the context of the detailed non-linear coagulation model of Fogelson et al., 10 a local analysis, thus the gradient information, proved sufficient for parameter ranking. 19 The gradient information is also useful when some of the inputs must be optimised by minimising the distance to some reference data. This path was for example followed in Fogelson et al. 10 to obtain realistic values of the kinetics of the sites where thrombin binds to activated platelets. In the same spirit, constant rates may require some adjustment by optimisation since the Michaelis-Menten kinetic often used to characterise enzymatic reactions may not be accurate when an enzyme or substrate is involved in more than one reaction. 19 Since the most efficient optimization algorithms are gradient-based (e.g., Newton-Raphson algorithm), it is clear that any efficient and accurate procedure to compute the gradient of the outcome would be useful.
The gradient of any outcome of a coagulation model can be computed by finite differences, which requires computing the outcome for many situations, at least n p þ 1 if is the total number of input parameters. On top of being time consuming when n p is getting large (e.g., n p ¼ 79 in Butenas et al. 7 ; n p ¼ 121 in Fogelson et al. 10 ), this approach has the disadvantage to rely on an increment, which is a user defined quantity used to perturb the input parameters one at a time. The coagulation model is then solved twice, once with the baseline value of the parameter of interest, once with its perturbed value, and the derivative of the output with respect to this parameter of interest is assessed from the output-to-parameter changes ratio. Since this ratio approximates the derivatives to first order in the increment, the temptation is to choose very small increment values. The risk is then that the numerical errors associated with each resolution of the coagulation model corrupts the change in the outcome, thus the derivative estimate. The best value of the increment in terms of numerical accuracy is often the result of a trial and error process and is problem specific. Thus, being able to faithfully assess the gradient of the outcome without defining an increment a priori would be advantageous.
Adjoint methods are efficient tools for analysing either the sensitivity or the receptivity of physical systems. 26 Sensitivity refers to the way the outcome of the system of interest depends on its input parameters whilst receptivity informs about how the system responds to external forcing. Using the adjoint formalism allows drastic computational reduction in situations where the number of system outputs is (much) smaller than the number of inputs. Note that this is the case for all the coagulation models where the number of inputs is easily of order 100 and the thrombin production is often abstracted by only a few quantities like the Endogeneous Thrombin Potential (ETP), maximum of thrombin generation, the clotting time and the maximum rate of thrombin production. 17 Adjoint methods have been used for dealing with different problems related to instabilities in fluid mechanics [27][28][29] but to the best of the author's knowledge, they have never been applied to coagulation models. The objective of this paper is to investigate if and how the adjoint formalism can be used to produce fast and accurate assessment of the gradient of the outcome of coagulation models with respect to the initial conditions and constant rates. The formalism is first established in Section 2 for a generic model and a pedagogical example is described in Section 3 for more clarity. The methodology is then applied to a detailed model for the extrinsic pathway in Section 4 to illustrate the efficiency of the approach.

| ADJOINT FORMALISM
Let us assume that the coagulation cascade is represented by the following set of first order differential equations: where Y is the column vector containing the molar fraction of the n s reactive species involved in the model and N k is a (non-linear) operator acting on Y. The index k reflects the fact that N k depends on the vector k which contains all the constant rates of the reactions forming the coagulation model. Once completed by the initial condition Y 0 ð Þ ¼ Y 0 , Equation 1 can be solved numerically to obtain the time evolution Y t ð Þ of each reactive species over any time range 0,T ½ . Now let us assume that the quantity (metric) J is used to characterise the outcome of Equation 1 (e.g., the total production of thrombin). Formally, this quantity depends on Y t ð Þ which itself depends on the kinetic constants k and initial conditions, globally referred to as p, a vector of size n p . The quantity of interest can then be seen as a function of Y t ð Þ and p, namely, J ¼ J Y,p ð Þ. Anyone interested in either the sensitivity of the model in Equation 1 or the optimal values of its parameters may want to assess the gradient of J with respect to p. From the chain rule, this gradient equals: In this expression, the partial derivatives of J are row vectors which can be easily assessed. However, the matrix dY=dp (which is made of n s rows and n p columns) is not known and is thus classically computed from finite differences. To avoid this costly and potentially inaccurate operation, it proves useful to introduce the following inner product: defined for any column vectors f and g of same length, as well as the following Lagrangian: where Y þ and λ are the yet unknown vectors of adjoint concentrations and Lagrange multipliers, respectively, and where the superscript t denotes the transpose of any matrix or vector. Obviously, ℒ is nothing but J for any solution of Equation 1 with initial condition Y 0 ð Þ ¼ Y 0 so that computing dJ =dp means in fact assessing dℒ=dp. Remarking that for any metric J which does not depend on time the following identity holds: the gradient of the Lagrangian can be expressed as follows, noting that the derivations with respect to p and time are independent operations: Integrating the d dt term by parts and gathering all the dY=dp then leads to: Inspecting Equation 6 shows that the adjoint variable and Lagrange multiplier can be selected in such a way that the terms involving dY dp cancel out in the previous equation. Indeed, if Y þ and λ are such that: the gradient of ℒ (and J ) is simply: Thus the gradient of the quantity of interest J with respect to all the initial conditions and constant rates can be computed from the knowledge of Y t ð Þ and Y þ t ð Þ only. The minus sign in the left hand side of Equation 7 indicates that the adjoint variable evolves backward in time, thus the initial condition at t ¼ T in Equation 7 and the use of the outcome Y þt 0 ð Þ in Equation 9. Note also that Equation 7 is linear in the adjoint by construction and that the ∂N k ∂Y term, namely the Jacobian of the direct problem Equation 1, may require to store the whole time evolution of Y t ð Þ and to interpolate in time during the resolution of the equation for the adjoint. At last, the quantity of interest mostly appears ( ∂J ∂p is most often zero in practical cases) through the forcing term ∂J ∂Y in Equation 7. This term is nothing but a row vector of length n s whose components are the functional derivatives of J with respect to the reactive species concentrations. It is current practice when dealing with coagulation to focus on the behaviour of thrombin, looking at the Endogeneous Thrombin Potential 3 (ETPused in the clinical routine as a marker of hyper-and hypo-coagulability), the instant where it starts increasing (referred to as the clotting time), its maximum rate and peak values. All the component of the vector ∂J ∂Y are then zero, except the one corresponding to thrombin, ∂J ∂y IIa , which may be assessed from the Fréchet derivative. The expressions of this functional derivative for a few classical J functionals are gathered in Table 1.

| A PEDAGOGICAL EXAMPLE
To clarify the formal developments made in Section 2, a simple case is considered in this section where the kinetic scheme is simply the combination of two species A and B which form the complex C: In the coagulation cascade, this type of reaction is for example relevant to the formation of prothrombinase from the activated factors X and V, or to the formation of tenase from activated VIII and IX.
The concentrations of the reactive species A, B and C are then solution of: where k is the constant rate of Equation 10. Noting y A0 , y B0 and y C0 the initial concentrations, the analytical solution of Equation 11 is then, assuming y A0 ¼ y B0 for simplicity: Two quantities of interest (metrics) will be considered in the following subsections.
The first quantity of interest is the total production of species C, namely, J ¼ R T 0 y C dt. From Equation 12, this quantity is then: and its derivatives read: Note the 1 2 term in the expression for ∂J ∂y A0 which comes from the fact that Equations 12 and 13 only hold for y A0 ¼ y B0 so that deriving the right hand side of Equation 13 with respect to y A0 actually gives ∂J In what follows, the adjoint framework detailed in Section 2 is used to retrieve the results in Equation 14. From Equation 7, the adjoint problem associated to Equations 11 and 13 reads: From Equation 15, the adjoint variable y þ C is simply Þwhilst after some algebra one may show that which indeed satisfies the "initial" condition y þ A T ð Þ ¼0 (recall that the adjoint variables evolve backward in time from T to 0). Assessing dJ dy A0 as y þ A 0 ð Þ=T in agreement with Equation 9 allows retrieving the first row of Equation 14. In the same way, since y þ C ¼ T TÀ t ð Þin the present case, one has y þ C 0 ð Þ ¼ T 2 , thus y þ C 0 ð Þ=T ¼ T, consistently with the y C0 derivative of Equation 13. Regarding the derivative of J with respect to k, Equations 9 and Equation 11 indicate that it equals: Injecting Equation 16 and y þ C ¼ T TÀ t ð Þinto Equation 17 gives, after some algebra, the second row of Equation 14.

| Time to threshold value
As a further illustration of the method, let us consider another quantity of interest, namely the time τ at which y C reaches a threshold value y th C . From Table 1, the T term in the last row of Equation 15 should be replaced by Þwhere τ is such that y C τ ð Þ ¼ y th C . From the third row of Equation 12, one has simply: and the equation for y þ C becomes From Equations 12 and 19 and the initial condition y þ C T ð Þ¼ 0, one may obtain that y þ C is zero in the interval τ, T ½ whilst it equals in the time window 0,τ ½ . Inserting this result into Equation 15 allows deriving the following expression for the other adjoint variables: The following results for the partial derivatives of τ follow: One may easily check that the above equations are retrieved by direct derivation of Equation 18. Note that the remark made in Section 3.1 about the extra 1 2 term that should be included because of the fact that Equations 12 and 13 only hold for y A0 ¼ y B0 remains valid.
Finally, the adjoint-based methodology allows computing the derivatives of any time independent quantity of interest (two were exemplified in this Section) with respect to the whole set of parameters of the simple kinetic scheme in Equation 10. One may also remark that the adjoint methodology gives access to exact values of the gradients, which is an advantage compared to the finite differences approximation.
The time evolution of species A and C are displayed in Figure 1 (y B is not shown since it equals y A ) together with the corresponding adjoint variables for the two quantities of interest considered. The analytical solutions derived above are in very good agreement with the numerical solutions (obtained from the Python 3 solve_ivp method) of Equations 11 and 15. Interpreting adjoint variables is not possible in general; in the present case one may remark that the fact that y þ A ¼ y þ C ¼ 0 for t > τ in the second case considered is coherent with the causality principle: changing either y A or y C after the instant where the threshold value of y C has been reached cannot modify the value of this instant.
The number of species and reaction rates is moderate in this particular case so that the gain in terms of computational burden is marginal. A more relevant and demanding case is thus considered in the next Section.

| APPLICATION TO A DETAILED COAGULATION MODEL
The adjoint-based methodology is now applied to a realistic kinetic scheme for the extrinsic pathway of the coagulation cascade. The one by Butenas et al., 7 slightly modified as in Danforth et al., 17,18 was selected for this purpose because it has already been used in parametric studies based on the Monte Carlo approach to assess how the generation of thrombin is affected by changes in the reactions rates 18 and the initial concentration of some of the coagulation factors. 17 Note however that the present adjoint formulation can be applied to any coagulation model which mathematically translates into a system of ordinary differential equations; this includes the schemes of Hockin et al., 6 Chatterjee et al., 8 Fogelson et al., 10 Susree and Anand, 24 Chen and Diamond, 15 amongst others.  Table 2.

| Baseline case
Solving the corresponding set of differential equations thanks to the Python 3 solve_ivp method (see Appendix A for more details), the global quantities relevant to the production of thrombin can be assessed; the corresponding values are gathered in Table 3 whilst Figure 2 displays the time evolution of the concentration of thrombin. They are in general in very good agreement with those reported in Danforth et al., 17 with less than 1% mismatch. Still, the discrepancy is larger for the clotting time, 4:4 mn = 264 s in Danforth et al. 17 against 237 s in the present study. One plausible explanation is a difference in the numerical methods at play. In the present case, the BDF stiff method as implemented in F I G U R E 1 Time evolutions of relevant quantities for the model Equation 10 obtained analytically (lines) and numerically (symbols). Left: Concentration of species A (solid) and C (dashed). Right: Adjoint variables A þ (solid) and C þ (dashed) for the metrics R T 0 y C dt (black) and τ such that y C τ ð Þ ¼ y th C (red). The results shown correspond to: T A B L E 2 Initial concentrations of the chemical species present at initial time. The initial concentration of any other species is zero. solve_ivp was first used to solve the coagulation scheme with 10 À10 relative and absolute error tolerances; then the result was interpolated thanks to a fifth order spline and the clotting time was deduced by solving y IIa ¼ 2 nM.

| Finite differences and adjoint-based gradients
To establish the validity of the adjoint-based gradients in the case of a complex coagulation scheme, a few derivatives of four metrics quantities relevant to thrombin are compared in Table 4. The finite difference approximation of the quantity J with respect to the parameter p i was obtained by adding either Δp ¼ 10%, Δp ¼ 1% or Δp ¼ 0:1% of the nominal value p i,0 of the parameter and forming the following growing ratio: T A B L E 3 Values of the global quantities relevant to thrombin. See Table 1 for the analytical expressions.  Table 2. The quantities of interest displayed in Table 3 are also illustrated in the Figure; the Endogeneous Thrombin Potential (ETP) is the area under the curve.
T A B L E 4 A selection of a few gradients of the metrics J 1 -J 4 defined in Table 3 with respect to arbitrarily selected initial concentrations and reaction rates. k 36 , k 1 , k 18 and k 2 are the reaction rates of TF = VIIa = Xa = TFPI ! TF = VIIa = Xa + TFPI, TF = VII ! TF + VII, IXa = VIIIa ! VIIIa + IXa and TF + VIIa ! TF = VIIa, respectively. The gradients estimates from finite differences (FD) have been computed as in Equation 24.
where J p 0 ð Þ is obtained by solving the coagulation model for the nominal values of the parameters (see Section 4.1), J p 0 þ Δp ð Þby adding Δp to the p i parameter only and solving the coagulation model once more. The values of the gradients displayed in Table 4 range over 18 decades and their units take 6 different values. Still, the adjoint-based assessment is consistent with the finite differences approximation in all the cases. Note that the discrepancy observed when Δp ¼ 10 % of p i,0 , and to a lesser extend when Δp ¼ 1 %, is due to the error associated with the first order accurate approximation in Equation 24 and virtually disappears when Δp ¼ 0:1 % of p i,0 .
T A B L E 5 Values of the scaled gradients of J 1 , J 2 , J 3 and J 4 with respect to the initial concentrations of the chemical species (left) and the reaction rates of the chemical reactions (right). For each metric J j and parameter p i , the scaled gradient dJ j dp i Ã is computed as in Equation 25. See Table 3 for the baseline values of the J 's and the present Table for the baseline parameters. Scaled gradients are displayed in light grey when the modulus in less than 0:2 and in bold when it is larger than 0:5. For each species and constant rate, the rank is based on the score S i obtained by averaging the absolute values of the four gradients as detailed in Equation 26. The chemical species not initially present in the baseline case are not ranked since they all correspond to p i,0 ¼ S i ¼ 0; for these species, the score values reported in the Table  were obtained  The complete lists of gradients of the four metrics J 1 , J 2 , J 3 , J 4 computed from the adjoint-based approach are displayed in Table 5. The derivatives with respect to the initial concentrations of the species are gathered in the left panel whilst the gradients with respect to the rate constants are shown on the right hand side part. Each entry features the corresponding gradient scaled by the baseline values of the metric of interest and of the parameter with respect to which the derivative is taken: Thus, a unitary entry means that a given percentage of change in the parameter induces the same percentage of change in the metric; conversely, small values correspond to cases where the metric is not very sensitive to the given parameter. This scaling allows comparing effects of parameters of different nature (modulus and/or units) but has the drawback to generate a zero entry for any species not present at the initial time in the baseline case (that is all the species except those in Table 2). Still, (some of) these species could potentially have an impact on the generation of thrombin. Thus, the scaling rule was changed for the species not initially present, a typical small value of 1 pM being used in place of their (zero) initial concentration.
The left hand side panel of Table 5 indicates that prothrombin (II) and anti-thrombin (ATIII) are by far the two species whose initial conditions impact the most the total production of thrombin (J 1 ). The same is true for the peak of thrombin metric, J 2 , if one considers only those species with non-zero initial concentration. In both cases, the same two species were identified in. 17 The tissue factor (TF) and tissue factor pathway inhibitor (TFPI) are, amongst the species initially present, the most prone to have an impact on the clotting time (J 3 metric). Only TFPI was identified by Danforth et al. 17 where TF was hold constant. The present results are again coherent with the previous findings 17 when considering the maximum rate of production of thrombin (metric J 4 ), with II, ATIII and TFPI being the most active species. No comparison with the previous work is possible when dealing with the species initially not present since their effect were not assessed. 17 None of these species would play a key role regarding ETP (J 1 ), but the peak thrombin may be impacted by the prothrombinase complex Xa = Va (combined or not with prothrombin), as strongly as by II and ATIII. The situation is the same for J 4 which also responds strongly to variation of the activated complex TF = VIIa (combined or not with factors IV, X and Xa). At last, one may note the very strong sensitivity of clotting time (J 3 ) on any initial concentration of the tenase complex IXa = VIIIa, and to a lesser extend to traces of prothrombinase.
Endogeneous thrombin potential, the peak value as well as the maximum rate of thrombin are strongly impacted by the rate constant of the neutralisation of thrombin by ATIII (reaction 40 in the right panel of Table 5). Conversely, the sequestration of the TF = VIIa = Xa complex by TFPI (reaction 35) has a limited impact on the total production but affects the clotting time, on top of the peak value and the maximum rate of production. This is in agreement with the findings of Danforth et al. 18 who show (see Figure 2 in this paper) that modifying the action of ATIII mostly modify the amplitude of the thrombin time evolution, whilst changing the TFPI action has also an impact on the time scale of the thrombin response.

| Parameter ranking
In order to provide a global ranking of the different input parameters with respect to the effect that their uncertainties have on thrombin generation, Danforth et al. 17,18 analysed the time evolutions of IIa as computed in their Monte Carlo database. To do so, they first computed for each parameter (either initial concentration or constant rate) the time evolution of the standard deviation of the thrombin concentration when one single parameter takes different values over its normal range. Then, from the time averaged standard deviations (one for each parameter), they computed the part of the total variance of IIa that is being "explained" by each parameter. This process can identify the constant rates for which variation in their values has the greatest consequences on thrombin generation but requires the prior generation of a Monte Carlo database.
Link et al. 19 developed a finite differences-based OAT local sensitivity analysis and concluded that this approach was sufficient to obtain good results in terms of parameter ranking. This idea is further tested in what follows, as an illustration of the usefulness of the knowledge of the gradient. The approach tested here is necessarily simpler and not as general/accurate as the one followed in Danforth et al. 17,18 since it relies only on the knowledge of the baseline case and its adjoint instead of a whole Monte Carlo database. To obtain a measure of the impact of each parameter on the thrombin generation process, the information contained in the gradient of each metric J 1 to J 4 should be used. There are an infinite ways of doing so; in the present study, the ranking of each parameter p i is generated from the score S i obtained by averaging the absolute value of the scaled gradients: The rank of each parameter is provided in the last column of the right (constant rates) and left (initial species concentrations) panels of ( Table 5). Although this procedure is very different from the Monte Carlo-based one, 17,18 they lead to very similar trends. For example, the 5 most influential constant rates in the right panel (k 40 , k 11 , k 35 , k 2 , k 0 ) are actually the same as in, 18 although ordered slightly differently (k 40 , k 2 , k 11 , k 35 , k 0 ); looking further in the ranking, 8 out of the 10 first ranked constant rates in 18 also appear in the top 10 values in (Table 5). Regarding the species (note that only the species initially present were ranked in the left panel of Table 5), II, ATIII and TFPI are identified in both cases as species for which variation in their values has the greatest consequences. Factor VIII was ranked #3 from the Monte Carlo approach, 18 but only #8 according to Table 5. Note however that VIIa and TF were not ranked in Danforth et al. 18 and that the score of factor VII, ranked #7, is very close to the one of factor VIII (see Table 5). Nevertheless, this factor seems to be less influential according to the present study, which is most likely due to the differences between the Monte Carlo approach and the local analysis presented here. As expected, TF has also a strong influence on thrombin generation (its effect was not assessed previously 17 ). Finally, both analyses indicate that uncertainties in factors V, VII and X have smaller impact on thrombin generation.
Two comments may be put forward to mitigate the good agreement between the current local analysis and the global Monte Carlo approach. First, given the local nature of the gradient-based score given in Equation 26, the fact that similar trends in the parameters ranking have been found for this particular scheme and operating point does not mean that the same would hold true for another nominal point of the same model, or for another coagulation scheme. Second, since the Monte Carlo approach accounts for the interactions between the parameters, the ranking given in 18 could change (significantly or not) if VIIa and TF would also have been considered in this study.
Each dimensional partial derivative discussed so far gives access to the variation ΔJ of any metric J when one single parameter changes by the amount Δp, the other keeping their baseline value: From the definition of the differential of J , the total variation of this metric when all the input parameters change at a time can be assessed by combining the partial derivatives and changes of all the parameters: The summation can also be taken over the initial concentrations of the species only, or the constant rates only, to produce ΔJ spec and ΔJ rate respectively. Assuming that the i th parameter of baseline value p i,0 experiences a change T A B L E 6 Sum of the absolute values of the scaled derivatives in Table 5. For each metric, the left entry corresponds to the summation over the species with non-zero initial concentration, the right one is for the whole set of constant rates, the value in parenthesis corresponding to the summation over the 10 , AEX % of its nominal value), the maximun change of J is obtained by fixing the sign of Δp i equal to the one of dJ dp i . At the end, the maximum relative change of J is approximated as: where dJ dp i Ã is the i th partial derivative scaled as in Equation 25. The values of the sum of the j dJ dp i Ã j terms are reported in Table 6 for the four metrics J 1 À J 4 , the summation being performed over either the 10 non-zero species or the whole set of rate constants. The corresponding evolutions of the maximum relative change of each metrics are displayed in Figure 3 which also shows some data obtained by solving the biochemical problem with appropriate values of the input parameters.
The maximum relative change in the maximum rate of growth of thrombin (J 4 ) is almost 10 times (9.79) larger than the relative change of each constant rate. This value decreases to approx. 6 when only the 10 most effective constant rates are considered; the amplification (6.18) is then comparable to the one observed for the species (6.43). The same trends are observed for the clotting time (J 3 ), although starting from a smaller amplification (4.09 instead of 9.79). The relative change of the peak value of thrombin (J 2 ) is also significantly (5.29) amplified compared to the relative changes in the reaction rates; it is also dramatically amplified by the 10 species initially present (4.48). ETP, even if a global quantity, can experience significant relative changes as well. Note also that contrary to the other metrics, J 1 is F I G U R E 3 Maximun variation of the four metrics J 1 À J 4 obtained by applying the same relative difference to either all the species with non-zero initial concentration (black) or the whole set of constant rates (red). The lines correspond to the linear variation from Equation 28 whilst each symbol corresponds to a direct solving of the coagulation problem. potentially more impacted by changes in the species concentrations (3.53) than by changes in the constant rates (2.54). Figure 3 shows that the approximation in Equation 28 is very accurate in all the cases when the input parameter varies by AE5 %, even though the relative change of some metrics can be quite large over this interval (approx. AE50 % for J 4 ). For some cases, the linear approximation remains accurate for larger variations of the input parameters, e.g. up to AE30% variation of the constant rates as far as J 1 is concerned. In general, the direct computations show that the maximum metrics variations are non linear with respect to the input changes for variations larger than 5% in amplitude. Note that the relative changes of the J 's cannot decrease below the À100% since all the metrics are positive by construction. The convex shape of the different curves in the negative range was thus expected. However, the fact that the linear approximation Equation 28 always underestimate the actual relative changes in the positive range could hardly be anticipated.
Of course, a one-at-a-time linear sensitivity analysis based on gradients computed at a single nominal point cannot represent the non-linearity of complex systems and the fact that some discrepancies arise for large enough amplitude is expected. At the same time, the fact that Equation 28 offers a faithful representation of the direct solving of the coagulation model at small amplitudes indicates that the adjoint-based method proposed in the paper leads to accurate gradients of the metrics J 1 -J 4 .

| CONCLUSION
An efficient and accurate method to compute the gradients of the outcome (metric) of coagulation models was successfully developed, validated and tested. To do so, an adjoint set of linear differential equations is introduced, of the same size as the non-linear set of differential equations governing the time evolution of the biochemical species. Once computed, the corresponding adjoint species can be used to explicitly compute the derivatives of the metric of interest. The nature of the latter appears as a forcing term in the adjoint equations and modifies also the way the derivatives are constructed from the adjoint variables.
The method was deployed for four metrics often used to characterise the efficiency of the coagulation cascade in closed systems: the Endogeneous Thrombin Potential, the thrombin peak value, the clotting time and the maximum rate of growth of thrombin during the burst. The application to a detailed coagulation scheme showed that the derivatives obtained from the adjoint-based methodology are fully consistent with those computed by finite differences. The advantage of the method is twofold: 1. the computed derivatives are exact (assuming the differential equations are properly solved) whilst they are only first order accurate with respect to the arbitrarily selected increment when finite differences are used, 2. the number of sets of differential equations to solve is only n J þ 1, where n J is the number of metrics for which the derivatives are sought for, whilst it is at least n p þ 1 for the finite differences approach.
In the example considered in this paper, where the coagulation scheme has 79 input parameters and there are 4 metrics of interest, the finite differences approach requires 80 simulations, namely, 16 times more than the requirement for the adjoint-based approach. The reduction of the computational load is thus dramatic. Note that in the case where no spatial heterogeneities are considered, modelling coagulation means solving a set of ordinary differential equations, a task that can be done in a few seconds by using appropriate stiff solver. Thus, making use of the brute force for assessing sensitivities is tractable, since performing hundreds of simulations can be achieved in a few minutes only. The situation is different if spatial heterogeneities are present since the model then relies on partial differential equations which may require minutes/hours to be solved properly. In this case, using an adjoint-based approach could become the only option. By demonstrating how to build a proper adjoint formulation for the spatially homogeneous biochemical part, the present study paves the way towards more complete formulations where spatial gradients would be present.
The maximum relative variation of the four metrics of interest were assessed from the computed gradients, showing that the linear approximation holds till each input parameter varies by less than approx. 5%. For larger variations of the inputs (either the initial concentrations or the rate constants), non-linear effects are not negligible, suggesting that a second order adjoint-based analysis would be a useful extension of the present work.