\section{Introduction}\label{sec:Intro}
Thrombosis, which is defined as excessive formation of blood clot or thrombus, is a common pathology in several cardiovascular diseases \cite{Nagareddy2013}, and blood-wetted medical devices \cite{DalSasso2019,Jaffer2019,Labarrere2020}. 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 \cite{Colman2006}. 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 (\cite{Leiderman2011a,Menichini2016,Yazdani:2017,Xu2008,Kim2013,Xu2017,Yazdani2018}). 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 \cite{Sorensen1999,Storti2014,Danes2019}; hemodynamics that regulate the transport of biochemical species \cite{Storti2014b,Fogelson2015}; and/or coagulation reactions that lead to the formation of fibrin \cite{Anand2006b,Yesudasan2019,Hansen2019}.
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 \cite{Belyaev2018}. These uncertainties are related to the fundamental structure of the model \cite{}, 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 \cite{Nikishova2018,Sankaran2016,Steinman2018}. UQ analyses in thrombosis models have been performed to identify the sensitivity of input parameters \cite{Danforth2009,MendezRojano2019a,Link2020,Melito2020}. In most sensitivity studies, the global Morris method \cite{Campolongo2007} or the Sobol' indices method \cite{Saltelli2010} 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. \cite{Melito2020} performed a variance based sensitivity analysis on the thrombosis model of Menichini et al. \cite{Menichini2016} 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. \cite{Link2020b} 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. \cite{Link2020b}  the synthetic databases built to perform sensitivity analysis can serve as hypothesis generation tool for \textit{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 multi-constituent thrombosis model introduced by Wu et al. \cite{Wu2017}. A widely used platelet function assay PFA-100\textregistered~ (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\textregistered ~closure time distribution in three scenarios: normal platelet count, thrombocytopenia, and vWF deficiency. The predicted closure times distribution align well with clinical data.

\section{Methods \label{sec:Methods}}
\subsection{Baseline Thrombosis Model}
The thrombosis model used in the current work is a variation of the approach introduced by Wu et al. \cite{Wu2017} with the addition of vWF activity. Figure~\ref{Fig:ThrombosisDiagram} A shows how platelets, the central player of the model, deposit to artificial surfaces. Figure~\ref{Fig:ThrombosisDiagram} 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) \textbf{AP} Activated platelets, 2) \textbf{RP} Resting Platelets, 3) \textbf{AP$_d$} Deposited activated platelets, 4) \textbf{RP$_d$} Deposited Resting Platelets 5) \textbf{ADP} Adenosine Diphosphate, 6) \textbf{TxA$_2$} Thromboxane A2, 7)\textbf{TB} Thrombin , 8) \textbf{PT} Prothrombin, 9) \textbf{AT} Anti-Thrombin, 10) \textbf{\vWFc} vWF Coiled, and 11) \textbf{\vWFs} vWF stretched.

\begin{figure}[h!]
\centering
\includegraphics[width=0.8\textwidth]{Fig1.png}
\caption{\label{Fig:ThrombosisDiagram} A) Diagram of platelet activation and deposition in the thrombosis model. $k_{rpd,b}$ and $k_{apd,b}$ 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$_d$ 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 \vWFs amplifies the deposition rate of free-flowing platelets by $\alpha_{k}$ and increases the resistance of deposited platelets to shear cleaning by $\alpha_{\tau}$.}
\end{figure}

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

\begin{equation}\label{eq:CDR}
\frac{\partial c_i}{\partial t} =  \nabla \cdot (D_i \nabla c_i) - \boldsymbol{v}_f  \cdot \nabla c_i + r_i ,
\end{equation}

where: $c_i $ is the concentration of species $i$, $D_i$ is the diffusion coefficient, $\boldsymbol{v}_f$ is the velocity vector field, and $r_i$ is the reaction source term that accounts for biochemical interactions such as platelet activation by thrombin or TxA$_2$.

\subsubsection{vWF activity}
To account for the role of vWF in high shear stress thrombosis the original model of Wu et al.\cite{Wu2017} was amended to include the enhanced platelet aggregation and adhesion driven by vWF activity. The right-hand side of Figure~\ref{Fig:ThrombosisDiagram} shows the two inter-convertible states of vWF in the model: collapsed \vWFc and stretched \vWFs. In all simulations, only \vWFc is introduced at the inlet of the domain and \vWFs is produced by means of \vWFc unfolding. The presence of stretched \vWFs polymers has a local thrombogenic effect via two mechanisms: it amplifies the deposition rate of free-flowing platelets by scaling it by $\alpha_{k}$ (Eq~(\ref{eq:kamp})) and increases the resistance of deposited platelets to shear cleaning by scaling the shear-cleaning parameters, $\tau_{emb}$ and $\tau_{emb,b}$, by $\alpha_{\tau}$ according to Eq~(\ref{eq:tauAmp})

\begin{equation}
\label{eq:kamp}
	\alpha_{k}=\begin{cases}
		1, & \text{if \vWFs $\leq$ \vWFsCrit}\\
		\frac{\textrm{\vWFs}}{\textrm{\vWFsCrit}}, & \text{if \vWFs $>$ \vWFsCrit}
	  \end{cases}
\end{equation}

\begin{equation}
\label{eq:tauAmp}
	\alpha_{\tau}=\begin{cases}
		1, & \text{if \vWFs $\leq$ \vWFsCrit}\\
		\frac{D_e(\dot{\gamma})}{D^{l}_{e}}, & \text{if \vWFs $>$ \vWFsCrit},
	  \end{cases}
\end{equation}

where \vWFsCrit is set as 5\% of the total vWF concentration; $D_e$ is the Morse potential's well depth proposed by Yazdani et al. \cite{Yazdani:2017} that correlates the platelet adhesion force to the local shear rate, $\dot{\gamma}$, according to Eq~(\ref{eq:DeYazdani});

\begin{equation}
\label{eq:DeYazdani}
	D_{e}(\dot{\gamma})=D_{e}^{h}\left[\tanh\left(\frac{\dot{\gamma}-\dot{\gamma}_{vW\!F\!}}{1000}\right)+\frac{D_{e}^{l}}{D_{e}^{h}}+1\right],
\end{equation}

$D^{l}_{e}$ and $D^{h}_{e}$ determine the adhesive forces at low and high shear rates, respectively, and $\dot{\gamma}_{vW\!F\!}$ = 5500 s\textsuperscript{-1} is the critical shear rate for vWF unfolding and marks the transition from low-shear to high-shear regime in this expression.

\vWFc 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. \cite{Babcock2003} where they define the flow-type $\lambda$ as:

\begin{equation}
\label{eq:flowType}
	\lambda=\frac{\|\boldsymbol{S}\|-\|\boldsymbol{\varOmega}\|}{\|\boldsymbol{S}\|+\|\boldsymbol{\varOmega}\|}
\end{equation}

where  $\|\boldsymbol{S}\|$ and $\|\boldsymbol{\varOmega}\|$ 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, $\lambda=1$ corresponds to purely extensional flow, $\lambda=0$ to simple shear flow, and $\lambda=-1$ to solid-body rotation. vWF will experience strong unfolding in flows with $\lambda\geq\lambda_{crit}$, simple-shear tumbling in $-\lambda_{crit}\leq\lambda<\lambda_{crit}$, and remain collapsed if $\lambda<-\lambda_{crit}$.

In simple shear and near-shear flows, $-\lambda_{crit}\leq\lambda<\lambda_{crit}$, the unfolding rate, $k_{c\text{-}s}$ (collapsed-to-stretched conversion), follows a function proposed by Lippok et al. \cite{Lippok2016} 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.~(\ref{eq:kcsShear})

\begin{equation}
\label{eq:kcsShear}
	k^{shear}_{c\text{-}s}(\dot{\gamma})=\frac{k^{\prime}}{1+\exp\left(-\frac{\dot{\gamma}-\dot{\gamma}_{vW\!F\!}}{\Delta\dot{\gamma}}\right)}
\end{equation}

\begin{equation}
\label{eq:kscShear}
	k^{shear}_{s\text{-}c}=k^{\prime}
\end{equation}

where $\dot{\gamma}_{vW\!F\!}$ = 5500 s\textsuperscript{-1} is the critical (half-maximum) shear rate for vWF unfolding and $\Delta\dot{\gamma}$ = 1271 s\textsuperscript{-1} denotes the width of the transition. $k^{\prime}$ is the nominal state conversion rate, given by $k^{\prime}=1/t_{vWF}$, where $t_{vWF}$ = 50 ms is the vWF unfolding time. This timescale for vWF unfolding was reported in experiments by Fu et al.\cite{Fu2017a} and Brownian dynamics simulations by Dong et al. \cite{Dong2019,Kania2021} To reflect the tumbling behavior of the vWF in this regime, the stretched-to-collapsed conversion rate, $k_{s\text{-}c}$ is set to the nominal value $k^{\prime}$ (Eq.~(\ref{eq:kscShear})).

In extensionally-dominated flows with $\lambda\geq\lambda_{crit}$, strong unfolding occurs if the modified Weissenberg number, $Wi_{ef\!f\!}$, given by Eq.~(\ref{eq:WiEff}), exceeds the critical value $Wi_{ef\!f\!,crit}$. Then, the unfolding rate $k_{c\text{-}s}^{steady}$ scales with $Wi_{ef\!f\!}$ according to Eq.~(\ref{eq:kcsSteady}), while $k_{s\text{-}c}^{steady}$ is zero.

If $Wi_{ef\!f\!}$ falls below $Wi_{ef\!f\!,crit}$, the stretched \vWFs exhibits hysteresis and does not collapse back to \vWFc until $Wi_{ef\!f\!}<Wi_{ef\!f\!,hyst}$, as shown in Eq.~(\ref{eq:kscSteady}).

\begin{equation}
\label{eq:WiEff}
	Wi_{ef\!f\!}=\sqrt{\lambda}\left(\|\boldsymbol{S}\|+\|\boldsymbol{\varOmega}\|\right)\tau_{rel}
\end{equation}

\begin{equation}
\label{eq:kcsSteady}
	k_{c\text{-}s}^{steady}=\begin{cases}
		k^{\prime}\frac{Wi_{ef\!f\!}}{Wi_{ef\!f\!,crit}}, & \text{if $Wi_{ef\!f\!} \geq Wi_{ef\!f\!,crit}$}\\
		k^{shear}_{c\text{-}s}, & \text{if $Wi_{ef\!f\!} < Wi_{ef\!f\!,crit}$}
	  \end{cases}
\end{equation}

\begin{equation}
\label{eq:kscSteady}
	k_{s\text{-}c}^{steady}=\begin{cases}
		0, & \text{if $Wi_{ef\!f\!,hyst} \leq Wi_{ef\!f\!} < Wi_{ef\!f\!,crit}$}.\\
		k^{\prime}, & \text{if $Wi_{ef\!f\!} < Wi_{ef\!f\!,hyst}$}.
	  \end{cases}
\end{equation}

where $\tau_{rel}$ 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.\cite{Sing2011b} Assuming the monomer diffusion time of $1.02\times10^{-3}$ s for $\tau_{rel}$ in Eq.~(\ref{eq:WiEff}), $Wi_{ef\!f,crit}=0.316$ and $Wi_{ef\!f,hyst}=0.053$.

If $\lambda<-\lambda_{crit}$, \vWFc remains collapsed and \vWFs reverts to a globular state, so $k_{c\text{-}s}=0$ and $k_{s\text{-}c}=k^{\prime}$.

A detailed validation of the vWF-thrombosis model is provided in Zhussupbekov et al. \cite{Zhussupbekov2021} for multiple academic simple flow configurations and complex in-vitro clotting scenarios.

\subsubsection{Flow dynamics}
The pressure and velocity fields $p$ and $\boldsymbol{v}_f$ are obtained solving the Navier-Stokes equations considering an incompressible Newtonian fluid:

\begin{equation}\label{ec:mass}
\frac{\partial \rho_f}{\partial t} + \nabla \cdot (\rho_f \textbf{v}_f) = 0
\end{equation}

\begin{equation}\label{ec:momentum}
\rho_f \frac{D \textbf{v}_f}{D t} = \nabla \cdot  \textbf{T}_f + \rho_f \textbf{b}_f + C_2 f(\phi) (\textbf{v}_f - \textbf{v}_T)
\end{equation}

where $\textbf{T}_f$ is the stress tensor of the fluid described as:

\begin{equation}\label{ec:stress tensor}
\textbf{T}_f = [-p(1-\phi)]\textbf{I} + 2\mu_f (1-\phi)\textbf{D}_f
\end{equation}

where $\mu_f$ is the asymptotic dynamic viscosity of blood. The scalar field $\phi$ is introduced to represent the volume fraction of deposited platelets (thrombus). The fluid phase $\rho_f$ is defined as

\begin{equation}\label{ec:density}
\rho_f = (1-\phi)\rho_{f0}
\end{equation}

where $\rho_{f0}$ is the fluid phase, $\textbf{b}_f$ is the body force, $\textbf{v}_f$ is and $ \textbf{v}_T $ are the velocity of the fluid and thrombus phases, respectively. $C_2$ is the resistance constant assuming that deposited platelets are densely compact spherical particles (2.78 $\mu$m), $f(\phi)$ is the hindrance function. For a full description of the original thrombosis model the reader is referred to the work of Wu et al.\cite{Wu2017} and Sorensen et al.\cite{Sorensen1999}.

\subsection{Case of Study PFA-100\textregistered}
The platelet function analyzer PFA-100\textregistered ~(Siemens Erlangen, Germany) is a coagulation testing device that is used to assess the primary hemostasis response \cite{Harrison2002}. Figure~\ref{fig:PFA-100} shows a cross section sketch of the PFA-100\textregistered ~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 $\Delta p$ of 4000 Pa. Shear rates inside the orifice reach values up to 6000 s$^{-1}$ promoting vWF unfolding leading to preferential thrombus formation\cite{Mehrabadi2016} 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 \cite{Harrison2002}. 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 \cite{Harrison2002}). In addition, CT times are affected by a variety of hematological and pharmacological factors such as hematocrit and aspirin.
\begin{figure}[h!]
\centering
\includegraphics[width=0.28\textwidth]{Fig2.png}
\caption{\label{fig:PFA-100} Cross section of PFA-100\textregistered ~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.}
\end{figure}

\subsection{Baseline Thrombosis Simulation}\label{subsec:baselineThrombosis}
\begin{figure}[h!]
\centering
\includegraphics[width=0.65\textwidth]{Fig3.png}
\caption{\label{fig:PFA-SetUp} A) Numerical setup of PFA-100\textregistered ~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.}
\end{figure}

An axisymmetric thrombosis simulation of the PFA-100\textregistered ~assay is the baseline case for the uncertainty quantification analysis. Figure~\ref{fig:PFA-SetUp}A shows the setup of the simulation. A uniform velocity profile $\textbf{v}_f = (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~\ref{tab:initial_conditions} 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 \cite{Steinlechner2009}. 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 where platelets can deposit. In addition, a diffusive flux of $ADP= 1 \times 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. \cite{Steinlechner2009}. All remaining variables were derived from the literature. The diffusion coefficient for each biochemical species were taken from \cite{Goodman2005}. We employed a Newtonian constitutive model for blood with a viscosity of 3.5 cP and a density of 1050 kg m$^{-3}$.

\begin{table}[h!]
\centering
\begin{tabular}{|c|c|}
\hline
Species & Concentration \\
\hline
RP &  216 $\times 10^3 $ Plt $\mu L^{-1}$ \\
AP &  2.16 $\times 10^3 $ Plt $\mu L^{-1}$ \\
\vWFc & 1000 nmol m$^{-3}$ \\
PT  &   1.1$\times 10^6 $ nmol m$^{-3}$  \\
TB  & 0 nmol m$^{-3}$ \\
AT  &  2.844  $\times 10^6 $ nmol m$^{-3}$ \\
ADP &  0 nmol m$^{-3}$ \\
TxA$_2$ &  0 nmol m$^{-3}$ \\
\hline
\end{tabular}
\caption{\label{tab:initial_conditions} Blood constituents inlet concentrations, the platelet count is representative of patients treated with ventricular assist devices \cite{Steinlechner2009}.}
\end{table}

\begin{figure}[h!]
\centering
\includegraphics[width=0.95\textwidth]{Fig4.png}
\caption{\label{fig:PFA-Observables} A) 3D rotational extrusion view of thrombus formation at 100 seconds of simulation, clot forms inside the membrane orifice. B) Quantities of interest ~($Y(X_i)$) for the UQ analysis, occlusion percentage across the top edge of the membrane (dashed line) and total thrombus volume inside (solid line).}
\end{figure}

The mesh was composed of 68650 hexahedral cells with a finest resolution of  $\Delta x = 3$ $\mu$m located in the membrane orifice. A time step of  $\Delta t = 1 \times 10^{-3} $ s was used to solve the flow and species conservation equations with an open source finite-volume software package (OpenFoam). 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~\ref{fig:PFA-SetUp}B and C.)

Figure~\ref{fig:PFA-Observables}A 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) \cite{Steinlechner2009}, the simulation closure time is defined at 80 \% of occlusion since full thrombotic closure is not numerically 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~\ref{fig:PFA-Observables}B.

\subsection{Uncertainty Quantification}
\subsubsection{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 and underlying numerical algorithmic and we will focus only on 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, we have selected a subset of 15 parameters that based on our experience have a significant impact on the thrombosis formation dynamics. The parameters are listed in Table~\ref{tab:parameters} 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~\ref{tab:parameters}.
\begin{table}[h!]
\centering
\begin{tabular}{| l | c | c |}
\hline
Parameter $X_i$ & Definition & Distribution (mean, variance)\\
\hline
1) RP   & Resting platelets & $\mathcal{N} (216e3,21.6e3)$ Plt $\mu L^{-1}$ \\
2) AP & Activated platelets & $\mathcal{N}(2.16e3,0.216e3)$ Plt $\mu L^{-1}$ \\
3) \vWFc & vWF concentration & $\mathcal{N}(1e3,0.1e3)$ nmol m$^{-3}$ \\
4) $K_{rpd}$ & Deposition rate RP (to membrane) & $\mathcal{N} (1.0e\text{-}8, 1.0e\text{-}9)$ m s$^{-1}$ \\
5) $K_{apd}$ & Deposition rate AP (to membrane) & $\mathcal{N}(3.0e\text{-}6, 3e\text{-}7)$ m s$^{-1}$ \\
6) $K_{aa}$ & Deposition rate of AP to $AP_d$ & $\mathcal{N}(3.0e\text{-}5, 3e\text{-}6)$ m s$^{-1}$ \\
7) $K_{ra}$ & Deposition rate of RP to $AP_d$ & $\mathcal{N}(3.0e\text{-}6,3e\text{-}7)$ m s$^{-1}$ \\
8) $\dot{\gamma}_{vW\!F\!}$ & critical shear rate for vWF unfolding & $\mathcal{N}(5.5e3,0.55e3)$ s$^-1$ \\
9) $\Delta\dot{\gamma}$  & vWF shear unfolding transition width & $\mathcal{N}(1.271e3, 0.1271e3)$ s$^2$ \\
10) $\lambda_{crit}$ & Flowtype critical value for strong unfolding & $\mathcal{N}(5e\text{-}3, 5 e\text{-}4)$ dimensionless \\
11) $t_{vWF}$ & vWF relaxation time & $\mathcal{N}(0.05, 0.005)$ s \\
12) $Wi_{ef\!f\!,crit}$ & $Wi_{ef\!f\!}$ critical value for strong unfolding & $\mathcal{N}(0.316, 0.0316) $ dimensionless \\
13) $Wi_{ef\!f\!,hyst}$ & $Wi_{ef\!f\!}$ hysteresis value & $\mathcal{N}(0.053,0.0053)$ dimensionless \\
14) \vWFsCrit & Critical \vWFs concentration value & $\mathcal{N}(50, 5)$\\
15) $D^{h}_{e}/D^{l}_{e}$ & Ratio of adhesive forces at high and low shear rates & $\mathcal{N}(500,50)$ dimensionless \\
\hline
\end{tabular}
\caption{\label{tab:parameters}Fifteen (15) parameters of the thrombosis model.}
\end{table}

\subsubsection{Polynomial chaos surrogate modeling framework}
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 $Y$ as accurately as possible.  It is constructed from a finite set of the thrombosis model predictions for 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~\cite{Rasmussen2006}, support vector machines~\cite{Smola1998}, stochastic interpolation~\cite{Tatang97}, 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 \cite{ghanem2003stochastic,le2010spectral,xiu2002wiener}. Let $(\Omega, \mathcal{B}, \mathcal{P})$ be the probability space where $\Omega$ is the space of random events $\omega$, this domain has a $\sigma$-algebra $\mathcal{B}$ and is equipped with a probability measure $\mathcal{P}$. The vector of normally-distributed random parameters can be written as ${\bX} \equiv {\bX}(\omega)=(X_{1}, \ldots, X_{d=15})$. If we consider, at a given time instant, the $d-$variate and second-order random variable $Y: \mathcal{I}_{{\bX}} \subseteq \mathbb{R}^d \rightarrow \mathbb{R}$, then $Y(t,{\bX}) \in \ell_2(\Omega, \mathcal{B}, \mathcal{P})$, can be expressed as the following expansion \cite{xiu2002wiener}:
\begin{align}
    Y(t,{\bX}) = \sum_{j=0}^{\infty} {{y}}_j(t)\, \mathrm{\psi_j}({\bX}),
    \label{gpc1}
\end{align}

where $\mathrm{\psi_j}({\bX}) = \prod_{i=1}^{d} \mathit{\psi_j^{(i)}}(X_{i})$ are the multivariate basis functions that form a complete basis, orthonormal with respect to the probability measure $\rho_{{\bX}}$ of the random input, and $\mathit{\psi_j^{(i)}}$ are the univariate Hermite polynomial functions along the $i^{\text{th}}$ dimension.
In practice, PCE expansion must be truncated, and writes:\footnote{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 $\Lambda_p$ is an index set for multi-index $\boldsymbol{\gamma} = (\gamma_1, \dots, \gamma_d) \in \mathbb{N}^d_0 $, then $\mathbb{P}_{\Lambda_p} \equiv \mathrm{span}\{{\psi}_{\boldsymbol{\gamma}} \  |\  \mathbf{\gamma} \in \Lambda_p\}$.
}
\begin{align}
    Y(t,\boldsymbol{\xi}) \approx \sum_{\boldsymbol{\gamma}  \in \Lambda_p} {y}_{\boldsymbol{\gamma} }(t) \,\mathrm{\psi_{\boldsymbol{\gamma} }}({\bX}),
    \label{truncatedPCE}
\end{align}

where ${y}_{\boldsymbol{\gamma} }(t)$ are the evolving {\it modal} coefficients corresponding to the $\mathrm{\psi_{\boldsymbol{\gamma} }}$ basis.
We will restrict ourselves to tensor-product polynomial spaces $\mathbb{P}_{\Lambda_p}$, where $\Lambda_p$ is an index set of ``degree" $p$, and where $P = \mathrm{dim}({\mathbb{P}_{\Lambda_p}})\equiv \#\Lambda_p$, 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 $q-$norm type of polynomial spaces, allowing for an hyperbolic truncation of the approximation basis (therefore reducing the computational cost):
         $\mathbb{P}_{\Lambda_{p,q}}$ with  index set
            $\Lambda_{p,q} = \{\by \in \mathbb{P}_{\Lambda_p} : \| \by \|_q = \left ( \sum_{\boldsymbol{\gamma}} y_{\boldsymbol{\gamma}}^q \right )^{1/q} \leq p \}$ \cite{Shen_SIAM2010}.   By default, the hyperbolic truncation with $q=1$ reduces to the standard total-degree truncation scheme. Decreasing the value of $q$ {\it de facto} decreases the number of polynomials of high interaction order kept in the expansion.

\subsubsection{Adaptively constructing the surrogate by least-squares minimization with regularization} \label{lsq}

Based on a training set $\{\bX^{(i)},Y^{(i)} \}_{i=1,\dots,N}$ collected at each time instant, ordinary linear regression can be used to compute the unknown coefficients $\boldsymbol{y} \equiv {y}_{\boldsymbol{\gamma}}$, \cite{choi2004polynomial,berveiller2006stochastic} by minimizing the residuals $\boldsymbol{r} \equiv \bY - \boldsymbol{\psi}_{\Lambda_p}
    \boldsymbol{y}$ in the $\ell_2-$norm through an optimization problem:
\begin{equation}
    \by = \argmin_{\boldsymbol{\mathrm{y}} \in \mathbb{R}^{P}} {\| \bY - \boldsymbol{\psi}_{\Lambda_p}
    \boldsymbol{\mathrm{y}}\|}_{2}^2,
    \label{regSol}
\end{equation}

where $\boldsymbol{\psi}_{\Lambda_p}$ is the measurement matrix corresponding to the gPC expansion in the index set $\Lambda_p$. The solution is obtained in matrix form as:
\begin{align}
    \by = \left( \boldsymbol{\Psi}_{\Lambda_p}^{\sf{T}} \boldsymbol{\Psi}_{\Lambda_p} \right)^{-1}\boldsymbol{\Psi}_{\Lambda_p}^{\sf{T}} \boldsymbol{\mathrm{Y}},
    \label{solLSQ}
\end{align}

where $\boldsymbol{\mathrm{Y}}$ is a vector of observations of size $N \times 1$, $\boldsymbol{\mathrm{\Psi}}_{\Lambda_p}$  the {\it measurement matrix} of size $N \times P$ with $\Psi_{ij}=\mathrm{\psi_j}({\bX_i})$, and $\boldsymbol{{y}}$ the vector of coefficients of size $P\times 1$. Unfortunately, when the number of simulations $N<P$, the problem is not well posed under this formulation and one needs to add a regularization term in the form of:
\begin{align}
    \boldsymbol{{y}} = \argmin_{\boldsymbol{\mathrm{y}} \in \mathbb{R}^{P}} \frac{1}{2} { \| \boldsymbol{Y} - \boldsymbol{\mathrm{\psi}}_{\Lambda_p} \boldsymbol{\mathrm{y}} \|}_{2}^2 + \lambda {\|\boldsymbol{\mathrm{y}}\|}_{1},
    \label{lassoDef}
\end{align}

with ${\|\boldsymbol{{y}}\|}_{1}=\sum_{{\boldsymbol{\gamma}  \in \Lambda_p}} \left | {y}_{\boldsymbol{\gamma}} \right |$, 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) \cite{efron2004least}, which has been adapted to the context of PCE by Blatman and Sudret\cite{BLATMAN20112345}. 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 \cite{van2016robust}.
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 $N$ metamodels, each one created on a reduced experimental design $\bX^l$, i.e. the training set obtained by removing a point $(\bX^{(l)},Y^{(l)})_{l \in \{1,2,\ldots,N\}}$, and comparing its prediction on the excluded scenario $\bX^{(l)}$ with the real value obtained from the thrombosis simulation. These $N$ errors are then averaged into a single value which may be collected for several time instants and different approximation bases for instance.

\subsubsection{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 $Y=f(t,\bX)$ is unique and hierarchic. We have:
\begin{equation}
Y=f(t,\bX)=\sum_{\bs\, \subseteq \{1,2,\ldots N_d\}} f_{\bs} (t,\bX_{\bs}),
\label{eq:ANOVA}
\end{equation}
where ${\bs}$ is a set of integers such that $\bX_{\bs}=(X_{s_1},\ldots, X_{s_N})$, with $N=\textrm{card}({\bs})=|{\bs}|$ and $f_{\varnothing}=f_0=\mathbb{E}[Y]$. In this way, the variance of the response can be decomposed as \cite{ANOVA_Efron81}:
\begin{equation}
\mathbb{V}(Y) \equiv \sigma^2(t)=\sum_{{\bs}\, \subseteq \{1,2,\ldots N_d\}} \sigma_{\bs}^2(t), \quad \textrm{with} \quad \sigma_{\bs}^2(t) = \mathbb{V}\left ( \mathbb{E}\left [ Y | \bX_{\bs} \right ] \right )- \sum_{ \substack{\bu \subset {\bs} \\ \bu \neq {\bs} \\ \bu \neq \varnothing} } \sigma_u^2(t),
\end{equation}
where $\mathbb{V}$ and $\mathbb{E}$ are  the variance and the expectation operators, respectively.
The normalized Sobol' indices  $S_{\bs}$ are defined as \cite{Sobol_1993}:
\begin{equation}
\quad S_{\bs} \equiv  \frac{\sigma_{\bs}^2(t)}{\sigma^2(t)}  \quad \textrm{and}  \sum_{ \substack{  \bs\, \subseteq \{1,2,\ldots N_d\} \\ \bs \neq \varnothing}} S_{\bs}(t) = 1,
\label{eq:sobol_coeff1}
\end{equation}
which measure the sensitivity of the variance of $Y$ due to the interaction between the variables
$\bX_{\bs}$, without taking into account the effect of the variables $\bX_{\bu}$  (for $\bu \subset \bs$ and $\bu \neq \bs$).
In this work, we consider first- and total-order Sobol' indices. The first-order quantifies the effect of the single parameter alone $X_s$ on the output. The total index of input variable $X_s$, denoted $S^T_s$, is the sum of all the Sobol' indices involving this variable: $S^T_s=\sum _{\bs\, \subseteq \{1,2,\ldots N_d\}} S_{\bs}$.
This latter definition is not practical since it would result in computing each index separately. Instead, by denoting $S_{\tilde{s}}$ the sensitivity measure of all the variables {\it excluding} variable $X_s$, the total index can be rewritten as:  $S^T_s=1-S_{\tilde{s}}$.
Previous Sobol' indices can be evaluated by Monte-Carlo estimators, but require very large samples. Instead, they can be obtained very straightforwardly directly from the polynomial chaos modal coefficients \cite{UQdoc_14_104}.

\section{Results}
\subsection{Multi-constituent Thrombosis Simulations}
The evolution in time of closure percentage and  thrombus volume is shown in Fig.\ref{fig:Observables} for 120 simulations. Closure percentage  (Fig.\ref{fig:Observables} 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.~\ref{fig:Observables}B 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 $m^3 \times 10^{-3}$.

Anticipating on the results of the following settings, both plots in Fig.~\ref{fig:Observables} also show the surrogate mean, std-based confidence interval, and a box plot of the 25$^{th}$ and 75$^{th}$ percentiles. It can be observed that the surrogate model is able to correctly reproduce the trend observed in the multi-constituent thrombosis simulations (gray plots).

\begin{figure}[!h]
\centering
\includegraphics[width=1\textwidth]{Fig5.png}
\caption{\label{fig:Observables} 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$^{th}$ and 75$^{th}$ percentiles of the surrogate model predictions. A) Occlusion percentage evolution in time B) Thrombus volume in the computational domain in time.}
\end{figure}

\subsection{PCE surrogate model training and validation}
A preliminary study was performed to determine a reasonable number of simulations to build our training database, achieving a balance between predictive accuracy and computational cost. To ensure that the parameter space $\bX$ was uniformly explored, a Sobol' quasi Monte Carlo sequence was used to generate a DoE matrix. Then we computed the LOO error as a function of the training sample size using as QoI the closure time as define in Section~\ref{subsec:baselineThrombosis}. (See Figure \ref{fig:LOOvsSample}.) This was repeated for five surrogate models of increasing complexity (increasing polynomial order $p$ with truncation parameter, $q=1$). With exception of the linear model, 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 $p=2$ and 3. Further increasing the polynomial order was found to deteriorate accuracy, reflecting the limited sample size. This motivated the use of an adaptive approach that searches, for each time step, a combination of polynomial orders $p \in \{1,6\}$ and truncation parameters $q \in \{0.5:0.1:1\}$. This yielded lower error than fixed polynomial order and truncation parameter approximations. In conclusion, this pilot study  demonstrated that 120 simulations were sufficient when combined with an adaptive polynomial approach.

\begin{figure}[h!]
\centering
\includegraphics[width=0.45\textwidth]{Fig6.pdf}
\caption{\label{fig:LOOvsSample} Time-averaged leave-one-out cross-validation error (LOO) of the Thrombus Volume parametric model as a function of the training sample size and model complexity. Adaptive construction of the surrogate model is optimal with a reasonable budget of 120 full thrombosis simulations.}
\end{figure}

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. To cross-validate the model considering time evolution, 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. Five predictions using this methodology are shown in Figure~\ref{fig:thrombusVolumeCrossValidation}, 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 \cite{MARELLI201867}]. Producing a single prediction from each PCE surrogate is not reasonable 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 $p_{\text{opt}}$), depending on the complex nature of the response to the parameters randomness along time.

\begin{figure}[ht]
\centering
\includegraphics[width=1\textwidth]{Fig7.png}
\caption{\label{fig:thrombusVolumeCrossValidation} Cross validation: predictions from PCE surrogate model (solid line), thrombosis simulations (circles), the shaded area represents the confidence intervals associated with the meta-model (2 $\pm$ std).}
\end{figure}

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 $t=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 $d=15$ dimensions and we only used a moderate computational budget of $N=120$ simulations.

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

\subsection{Sobol' Indices}
The first order and total Sobol' indices for the 15 parameters computed using the closure percentage as QoI are shown in Fig.~\ref{fig:Sobol1stOcclusion}. 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 $50s<t<200s$ which is in line with the fact that most variations in the QoI are observed within this time range (see Fig.~\ref{fig:Observables} A). The physical significance of the parameters identified by the first Sobol' indices is related to three main factors:
\begin{itemize}
\item \textbf{Thrombus constitution.} Resting platelet count (parameter 1) is the main constituent of the thrombus, therefore, it is understandably influential in the closure percentage.
\item \textbf{vWF activity.} The parameters 3, 8, 11, 12, and 14 ($vWF_c$, $\dot{\gamma}_{vWF}$, $t_{vWF}$, $Wi^{crit}_{eff}$, and $vWF^{crit}_{s}$) 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.
\item \textbf{Platelet-platelet adhesion.} The parameter 6 ($k_{aa}$) is the rate of platelet to platelet adhesion which directly influences the thrombus growth rate.
\end{itemize}

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.~\ref{fig:Sobol1stOcclusion} B. Significant effects for 11 out of 15 parameters are observed indicating non-negligible 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.

\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{Fig8.png}
\caption{\label{fig:Sobol1stOcclusion} Sobol' indices evolution over time for the QoI closure percentage. A) First order indices multiplied by the output PCE variance. B) Total Sobol' indices.}
\end{figure}

Figure~\ref{fig:ThrombusVolume} 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.~\ref{fig:Observables}B). 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 ($\dot{\gamma}_{vWF}$) 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 \ref{fig:ThrombusVolume}B.) 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 ($\Delta\dot{\gamma}$ and $k_{ra}$) show a low but non-negligible impact, these parameters are related to direct vWF activity and platelet deposition respectively.
\begin{figure}[h!]
\centering
\includegraphics[width=1\textwidth]{Fig9}
\caption{\label{fig:ThrombusVolume} Sobol' indices evolution over time for the thrombus volume QoI. A) First order indices. B) shows the total Sobol' indices.}
\end{figure}

\subsection{Clinical Application: PFA-100\textregistered ~Closure Time}\label{sec:ResultsClinic}
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\textregistered ~as QoI as defined in Section~\ref{subsec:baselineThrombosis}. 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~\ref{tab:PCE_Clinical} were drawn using the Monte Carlo sampling technique described previously. The 13 supplementary model parameters that do not appear in Table~\ref{tab:PCE_Clinical} were held constant, as defined in Table~\ref{tab:parameters}.

\begin{table}[h!]
\centering
\begin{tabular}{|c|c|c|}
\hline
\textbf{Case} & \textbf{Parameter} & \textbf{Distribution} \\
\hline
\textbf{Normal} & RP Resting platelet count & $\mathcal{N} (300e3,30e3)$ Plt $ \mu L^{-1}$  \\
\textbf{vWFD} & \vWFc (70 \% of vWF concentration) & $\mathcal{N}(700,70)$ nmol m$^{-3}$ \\
\textbf{Thrombocytopenia} & RP (50\% of Normal Resting platelet count)  & $\mathcal{N} (150e3,15.0e3)$ Plt $ \mu L^{-1}$ \\
\hline
\end{tabular}
\caption{\label{tab:PCE_Clinical} Definition of random parameters for the three cases.}
\end{table}

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~\ref{fig:PCEClossureTimes} (solid lines) superimposed with clinical data derived from 77 patients, reported by Harrison et al. \cite{Harrison2002} demonstrating remarkable agreement; although the simulated closure time for vWFD and thrombocytopenia lag the clinical data slighltly.

\begin{figure}[h!]
\centering
\includegraphics[width=0.7\textwidth]{Fig10.png}
\caption{\label{fig:PCEClossureTimes} Distribution of simulated and clinically reported PFA-100\textregistered ~closure time for normal blood, vWF disease and thrombocytopenia. The histogram corresponds to clinical data from \cite{Harrison2002} using ADP cartridges. The blue line is the probability density function computed with one hundred thousand realizations of the PCE surrogate model.}
\end{figure}

\section{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 \cite{Peach:2014,WanAbNaim2018,Yazdani2021}. 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. \cite{Melito2020} in which the authors found that only 4 of the 9 parameters in the model of Menichini \cite{Menichini2016} 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 \cite{Wu2017}, 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, 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 asyptotic value of approximatley 90\%  due to extremely high shear rates as the orifice diameter approaches zero. Finally, the clinical database for normal blood was constructed with n=50 patients; for the vWFD and thrombocytopenia, only n=15 patients and n=12 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 = $280e3$ Plt $\mu$ 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.

\section{Acknowledgments}
This research was supported by NIH R01 HL089456 and NIH R01 HL086918.

\bibliographystyle{unsrt}


\section{Introduction}

Cancer diagnosis and treatment plans are guided by multiple streams of data acquired from several modalities, such as radiology scans, molecular profiling, histology slides, and clinical variables. Each characterizes unique aspects of tumor biology and, collectively, they help clinicians understand patient prognosis and assess therapeutic options. Advances in molecular profiling techniques have enabled the discovery of prognostic gene signatures, bringing precision medicine to the forefront of clinical practice \cite{el-deiry_current_2019}. More recently, computational techniques in the field of radiology have identified potential imaging-based phenotypes of treatment response and patient survival. Such approaches leverage large sets of explicitly designed image features (commonly known as radiomics \cite{gillies_radiomics:_2015}) or entail the novel discovery of image patterns by optimizing highly parameterized deep learning models such as convolutional neural networks (CNN) \cite{saba_present_2019} for prediction. Along similar lines, the digitization of histopathology slides has opened new avenues for tissue-based assays that can stratify patients by risk from H\&E slide images alone \cite{skrede_deep_2020}. Given the complementary nature of these various modalities in comprehensive clinical assessment, we hypothesize that their combination in a rigorous machine learning framework may predict patient outcomes more robustly than qualitative clinical assessment or unimodal strategies.

Glioma is an intuitive candidate for deep learning-based multimodal biomarkers owing to the presence of well-characterized prognostic information across modalities \cite{louis_2016_2016}, as well as its severity \cite{siegel_cancer_2017}.  Gliomas can be subdivided by their malignancy into histological grades II-IV \cite{louis_2016_2016}. Grades differ in their morphologic and molecular heterogeneity \cite{olar_using_2014}, which correspond to treatment resistance and short-term recurrence \cite{stupp_radiotherapy_2005,parker_intratumoral_2016}. Quantitative analysis of glioma \cite{bae_radiomic_2018} and its tumor habitat \cite{prasanna_radiomic_2017} on MRI has demonstrated strong prognostic potential, as well as complex interactions with genotype \cite{beig_radiogenomic-based_2020} and clinical variables \cite{beig_sexually_2021}.

Most deep multimodal prediction strategies to date have focused on the fusion of biopsy-based modalities \cite{chen_pathomic_2019,mobadersany_predicting_2018,cheerla_deep_2019}. For instance, previous work integrating molecular data with pathology analysis via CNN or graph convolutional neural networks (GCN) has shown that a deep, multimodal approach improves prognosis prediction in glioma patients  \cite{chen_pathomic_2019,mobadersany_predicting_2018}. Likewise, Cheerla et al. integrated histology, clinical, and sequencing data across cancer types by condensing each to a correlated prognostic feature representation \cite{cheerla_deep_2019}. Multimodal research involving radiology has been predominantly correlative in nature \cite{beig_sexually_2021,beig_radiogenomic-based_2020}. Some have explored late-stage fusion approaches combining feature-based representations from radiology with similar pathology \cite{vaidya_raptomics_2018} or genomic features \cite{subramanian_multimodal_2020} to predict recurrence. While promising, these strategies rely on hand-crafted feature sets and simple multimodal classifiers that likely limit their ability to learn complex prognostic interactions between modalities and realize the full additive benefit of integrating diverse clinical modalities.

To our knowledge, no study to date has combined radiology, pathology, and genomic data within a single deep learning framework for outcome prediction or patient stratification. Doing so requires overcoming several challenges. First, owing to the difficulty of assembling multimodal datasets with corresponding outcome data in large quantities, fusion schemes must be highly data efficient in learning complex multimodal interactions. Second, the presence of strongly correlated prognostic signals between modalities \cite{cheerla_deep_2019} can create redundancy and hinder model performance.

In this paper, we introduce a deep learning framework that combines radiologic, histologic, genomic, and clinical data into a fused prognostic risk score. Using a novel technique referred to as Deep Orthogonal Fusion (DOF), we train models using a Multimodal Orthogonalization (MMO) loss function to maximize the independent contribution of each data modality, effectively improving predictive performance. Our approach, depicted in Fig. \ref{fig1}, first trains unimodal embeddings for overall survival (OS) prediction through a Cox partial likelihood loss function. Next, these embeddings are combined through an attention-gated tensor fusion to capture all possible interactions between each data modality. Fusion models are trained simultaneously to predict OS and minimize the correlation between unimodal embeddings. We emphasize the following contributions:

\textbf{Deep Fusion of Radiology, Pathology, and Omics Data}: We present a powerful, data-efficient framework for combining oncologic data across modalities. Our approach enabled a previously unexplored deep integration of radiology with tissue-based modalities and clinical variables for patient risk stratification. This fusion model significantly improved upon unimodal deep learning models. In particular, we found that integrating radiology into deep multimodal models, which is under-explored in previous prognostic studies, conferred the single greatest performance increase. This finding suggests the presence of independent, complementary prognostic information between radiology and biopsy-based modalities and warrants their combination in future prognostic studies.

\textbf{MMO}: To mitigate the effect of inherent correlations between data modalities, we present an MMO loss function that penalizes correlation between unimodal embeddings and encourages each to provide independent prognostic information. We find that this training scheme, which we call DOF, improves prediction by learning and fusing disentangled, prognostic representations from each modality. DOF was also found to outperform a fusion scheme that enforces correlated representations between modalities \cite{cheerla_deep_2019}, emphasizing that the dissimilarity of these clinical data streams is crucial to their collective strength.

\textbf{Multi-parametric Radiology FeatureNet:} A neural network architecture that can fuse CNN-extracted deep features from local tumor regions on multiple image sequences (e.g., Gd-T1w and T2w-FLAIR scans) with global hand-crafted radiomics features extracted across the full 3D region-of-interest.

\textbf{Independent prognostic biomarker of OS in glioma patients:} Using 15-fold Monte Carlo cross-validation with a 20\% holdout test set, we evaluate deep fusion models to predict glioma prognosis. We compare this multimodal risk score with existing prognostic clinical subsets and biomarkers (grade, \textit{IDH} status) and investigate its prognostic value within these outcome-associated groups.

\begin{figure}
\centering
\includegraphics[width=1\textwidth]{fusion_framework_V6.png}
\caption{DOF model architecture and training.} \label{fig1}
\end{figure}

\section{Methodology}

Let $X$ be a training minibatch of data for $N$ patients, each containing $M$ modalities such that $X = [x_1, x_2, ... , x_M]$. For each modality $m$, $x_m$ includes data from for $N$ patients. $\Phi_m$ denotes a trainable unimodal network, which accepts $x_m$ and generates a deep embedding $h_m = \Phi_m(x_m) \in \mathbb{R}^{l_1 x N}$.

\subsection{Multimodal Fusion}
When $M>1$, we combine embeddings from each modality in a multimodal fusion network. For each $h_m$, an attention mechanism is applied to control its expressiveness  based on information from the other modalities. An additional fully connected layer results in $h_m^S$ of length $l_2$. Attention weights of length $l_2$ are obtained through a bilinear transformation of $h_m$ with all other embeddings (denoted as $H_{\cancel{m}}$), then applied to $h_m^S$ to yield the attention-gated embedding:
\begin{equation}
 h_m^* = a_m * h_m^S = \sigma (h_m^T * W_A * H_{\cancel{m}})* h_m^S.
\end{equation}

To capture all possible interactions between modalities, we combine attention-weighted embeddings through an outer product between modalities, known as tensor fusion \cite{zadeh_tensor_2017}. A value of 1 is also included in each vector, allowing for partial interactions between modalities and for the constituent unimodal embeddings to be retained. The output matrix
\begin{equation}
 F =  \begin{bmatrix} 1 \\
h_{1}^* \end{bmatrix} \otimes \begin{bmatrix} 1 \\
 h_{2}^* \end{bmatrix} \otimes ... \otimes \begin{bmatrix} 1 \\
h_{M}^* \end{bmatrix}
\end{equation}
is an $M$-dimensional hypercube of all multimodal interactions with sides of length $l_2 + 1$. Fig. \ref{fig1} depicts $F$ for the fusion of radiology, pathology, and genomic data. It contains subregions corresponding to unaltered unimodal embeddings, pairwise fusions between 2 modalities, and trilinear fusion between all three of the modalities. A final set of fully connected layers, denoted by $\Phi_F$, is applied to tensor fusion features for a final fused embedding $h_F = \Phi_F (F)$.

\subsection{MMO Loss}
To address the shortcoming of multimodal models converging to correlated predictors, we introduce MMO loss. Inspired by Orthogonal Low-rank Embedding \cite{lezama_ole_2017}, we stipulate that unimodal embeddings preceding fusion should be orthogonal. This criterion enforces that each modality introduced contributes unique information to outcome prediction, rather than relying on signal redundancy between modalities. Each $\Phi_m$ is updated through MMO loss to yield embeddings that better complement other modalities. Let $ H \in \mathbb{R}^{l_1 x M*N}$ be the set of embeddings from all modalities. MMO loss is computed as
\begin{equation}
L_{MMO} = \frac{1}{M*N} \sum_{m=1}^M max(1, ||h_m||_*) - ||H||_*
\end{equation}
where $ || \cdot ||_*$ denotes the matrix nuclear norm (i.e., the sum of the matrix singular values). This loss is the difference between the sum of nuclear norms per embedding and the nuclear norm of all embeddings combined. It penalizes the scenario where the variance of two modalities separately is decreased when combined and minimized when all unimodal embeddings are fully orthogonal. The per-modality norm is bounded to a minimum of 1 to prevent the collapse of embedding features to zero.

\subsection{Cox Partial Likelihood Loss}
The final layer of each network, parameterized by $\beta$, is a fully connected layer with a single unit. This output functions as a Cox proportional hazards model using the deep embedding from the previous layer, $h$, as its covariates. This final layer's output, $\theta$, is the log hazard ratio, which is used as a risk score. The log hazard ratio for patient $i$ is denoted as $ \theta_i = h_i^T * \beta $.

We define the negative log likelihood $L_{pl}$ as our cost function
\begin{equation}
L_{pl} = - \sum_{i:  E_i=1} \left ( \theta_i - log \sum_{j: t_i \geq  t_j} e^{\theta_j}    \right )
\end{equation}
where $t \in \mathbb{R}^{Nx1} $ indicates the time to date of last follow up. The event vector, $E\in\{0,1\}^{Nx1}$, equals 1 if  an event was observed (death) or 0 if a patient was censored (still alive) at time of last follow up. Each patient $i$ with an observed event is compared against all patients whose observation time was greater than or equal to $t_i$. Networks are trained using the final loss $L$ which is a linear combination of the two loss functions specified above
\begin{equation}
 L = L_{pl} + \gamma L_{MMO}
\end{equation}
where $\gamma$ is a scalar weighting the contribution of MMO loss relative to Cox partial likelihood loss. When training unimodal networks, $\gamma$ is always zero. Performance for various values of $\gamma$ are included in the Table \ref{tab:tabs4}.

\subsection{Modality-specific Networks for Outcome Prediction}
\textbf{Radiology: }A multiple-input CNN was designed to incorporate multiparametric MRI data and global lesion measurements, shown in Fig. \ref{figS1}. The backbone of the network is a VGG-19 CNN \cite{simonyan_very_2015} with batch normalization, substituting the final max pooling layer with a 4x4 adaptive average pooling. Two pre-trained \cite{deng_imagenet_2009} CNN branches  separately extract features from Gd-T1w and T2w-FLAIR images, which are then concatenated and passed through a fully connected layer. A third branch passes hand-crafted features (described in section 3) through a similar fully connected layer. Concatenated embeddings from all branches are fed to 2 additional fully connected layers. All fully connected layers preceding the final embedding layer have 128 units.

\textbf{Histology, genomic, and clinical data:} We reused the models proposed in \cite{chen_pathomic_2019} - a pre-trained VGG-19 CNN with pretrained convolutional layers for Histology and a Self-Normalizing Neural Network (SNN) for genomic data. We also use this SNN for analysis of clinical data, which was not explored in \cite{chen_pathomic_2019}.

  \begin{figure}
 \centering
\includegraphics[width=.9\textwidth]{sampling_strategy_V3.png}
\caption{Sampling multiple radiology \& pathology images for patient level risk scores.} \label{fig2}
\end{figure}

\section{Experimental Details}

\textbf{Radiology}: 176 patients (see patient selection in Fig. \ref{figS2}) with Gd-T1w and T2w-FLAIR scans from the TCGA-GBM \cite{scarpace_radiology_2016} and TCGA-LGG \cite{pedano_radiology_2016} studies were obtained from TCIA \cite{clark_cancer_2013} and annotated by 7 radiologists to delineate the enhancing lesion and edema region. Volumes were registered to the MNI-ICBM standardized brain atlas with 1 mm isotropic resolution, processed with N4 bias correction, and intensity normalized. 96x96x3 patches were generated from matching regions of Gd-T1w and T2w-FLAIR images within the enhancing lesion. For each patient, 4 samples were generated from four even quadrants of the tumor along the z-axis. Patch slice position was randomized in unimodal training and fixed to the middles of quadrants during inference and fusion network training. Nine features including size, shape, and intensity measures were extracted separately from Gd-T1w and T2w-FLAIR images, and summarized in three fashions for a total of 56 handcrafted features, listed in Table \ref{tab:tabs1}.

\textbf{Pathology and Genomics:} We obtained 1024×1024 normalized regions-of-interest (ROIs) and DNA sequencing data curated by \cite{mobadersany_predicting_2018}. Each patient had 1-3 ROIs from diagnostic H\&E slides, totaling 372 images. DNA data consisted of 80 features including mutational status and copy number variation (Table \ref{tab:tabs2}).

\textbf{Clinical information: }14 clinical features were included into an SNN for the prediction of prognosis. The feature set included demographic and treatment details, as well as subjective histological subtype (see Table \ref{tab:tabs3}).

\textbf{Implementation Details:} The embedding size for unimodal networks, $l_1$, was set to 32. Pre-fusion scaled embedding size, $l_2$, was 32 for $M$=2, 16 for $M$=3, and 8 for $M$=4. Post-fusion fully connected layers consisted of 128 units each. The final layer of each network had a single unit with sigmoid activation, but its outputs were rescaled between -3 and 3 to function as a prognostic risk score. Unimodal networks were trained for 50 epochs with linear learning rate decay, while multimodal networks were trained for 30 epochs with learning rate decay beginning at the 10th epoch. When training multimodal networks, the unimodal embedding layers were frozen for 5 epochs to train the fusion layers only, then unfrozen for joint training of embeddings and fusion layers.

\textbf{Statistical Analysis: }All models were trained via 15-fold Monte Carlo cross-validation with 20\% holdout using the patient-level splits  provided in \cite{mobadersany_predicting_2018}. The primary performance metric was the median observed concordance index (C-index) across folds, a global metric of prognostic model discriminant power. We evaluated all possible combinations of a patient’s data (see sampling strategy in Fig. \ref{fig2}) and used the 75th percentile of predicted risk score as their overall prediction. C-indexes of the best-performing unimodal model and the DOF multimodal model were compared with a Mann-Whitney $U$ test \cite{mann_test_1947}. Binary low/high-risk groups were derived from the risk scores, where a risk score $>$0 corresponded to high risk. For Kaplan-Meier (KM) curves and calculation of hazard ratio (HR), patient-level risk scores were pooled across validation folds.

\begin{table}[]
\centering
\caption{Median C-index of unimodal and fusion models with and without MMO loss. }
\label{tab:table1}
\begin{tabular}{@{}cccc@{}}
\toprule
\textbf{Group}  & \textbf{Model}    & \textbf{Cox Loss Only} & \textbf{With MMO Loss} \\ \midrule
Unimodal        & Rad               & 0.718 ± 0.064          & --                     \\
                & Path              & 0.715 ± 0.054          & --                     \\
                & Gen               & 0.716 ± 0.063          & --                     \\
                & Clin              & 0.702 ± 0.049          & --                     \\
Pairwise Fusion & Path+Gen          & 0.711 ± 0.055          & 0.752 ± 0.072          \\
                & Gen+Clin          & 0.702 ± 0.053          & 0.703 ± 0.052          \\
                & Rad+Gen           & 0.761 ± 0.071          & 0.766 ± 0.067          \\
                & Rad+Path          & 0.742 ± 0.067          & 0.752 ± 0.072          \\
                & Rad+Clin          & 0.746 ± 0.068          & 0.736 ± 0.067          \\
                & Path+Clin         & 0.696 ± 0.051          & 0.690 ± 0.043          \\
Triple Fusion   & Path+Gen+Clin     & 0.704 ± 0.059          & 0.720 ± 0.056          \\
                & Rad+Path+Clin     & 0.748 ± 0.067          & 0.741 ± 0.067          \\
                & Rad+Gen+Clin      & 0.754 ± 0.066          & 0.755 ± 0.067          \\
                & Rad+Path+Gen      & 0.764 ± 0.062          & \textbf{0.788 ± 0.067} \\
Full Fusion     & Rad+Path+Gen+Clin & \textbf{0.775 ± 0.061} & 0.785 ± 0.077          \\ \bottomrule
\end{tabular}
\end{table}

\section{Results and Discussion}

Genomic- and pathology-only model performance metrics are practically similar (Table \ref{tab:table1}). However, the CNN-only (C-index=0.687 ± 0.067) and feature-only (C-index=0.653 ± 0.057) configurations of the radiology model underperform relative to the aforementioned unimodal models. Combining the radiology CNN features with the handcrafted features results in the strongest unimodal model. In contrast, clinical features are the least prognostic unimodal model.

\begin{figure}
    \centering
    \def\columnwidth{\columnwidth}
    \input{miccai_fig_3_cameraready_tightcrop.pdf_tex}
    \caption{Stratification by (a) grade, (b) IDH mutation status, and (c) DOF risk groups.} \label{fig3}
\end{figure}

 \begin{figure}
    \centering
    \def\columnwidth{\columnwidth}
    \input{miccai_fig_4_cameraready_tightcrop_tighter.pdf_tex}
    \caption{DOF risk groups stratify patients by OS within (a,b) grade \& (c,d) IDH subsets.} \label{fig4}
\end{figure}

Deep fusion models integrating radiology outperform individual unimodal models, naive ensembles of unimodal models, as well as fusions of only clinical and/or biopsy-derived modalities.  The full fusion model (C-index=0.775 ± 0.061) achieves the best performance when trained with Cox loss \cite{ching_cox-nnet_2018} alone, second only to the Rad+Path+Gen model trained with MMO loss. Naive late fusion ensembles (i.e., averaging unimodal risk scores) exhibit inferior performance for Rad+Path+Gen with (C-index=0.735 ± 0.063) and without (C-index=0.739 ± 0.062) clinical features, confirming the benefits of deep fusion.

The addition of MMO loss when training these deep fusion models consistently improves their performance at five different weightings (Table \ref{tab:tabs4}), with best performance for both at $\gamma=.5$. When all fusion models are trained at this weighting, 8 of 11 improve in performance. DOF combining radiology, pathology, and genomic data predicts glioma survival best overall with a median C-index of 0.788 ± 0.067, a significant increase over the best unimodal model (p=0.023).

An ablation study was conducted to investigate the contributions of components of the fusion module (modality attention-gating and tensor fusion). We found that a configuration including both yields the best performance, but that strong results can also be achieved with a simplified fusion module (Table \ref{tab:tabs5}).

In Fig. \ref{fig3}, KM plots show that the stratification of patients by OS in risk groups derived from this model perform comparably to established clinical markers. In Fig. \ref{fig4}, risk groups further stratify OS within grade and \textit{IDH} status groups. In sum, these results suggest that the DOF model provides useful prognostic value beyond existing clinical subsets and/or individual biomarkers.

To further benchmark our approach, we implemented the fusion scheme of \cite{cheerla_deep_2019}, who combined pathology images, DNA, miRNA, and clinical data, which we further modified to also include radiology data. The network and learning approach is described in-depth in Table \ref{tab:tabs6}. In contrast to DOF, \cite{cheerla_deep_2019} instead seeks to maximize the correlation between modality embeddings prior to prediction. A model combining radiology, pathology, and genomic data achieved C-index=0.730 ± 0.05, while a model excluding the added radiology arm stratified patients by OS with C-index=0.715 ± 0.05.

\section{Conclusions}
We present DOF, a data efficient scheme for the novel fusion of radiology, histology, genomic, and clinical data for multimodal prognostic biomarker discovery. The integration of multi-dimensional data from biopsy-based modalities and radiology strongly boosts the ability to stratify glioma patients by OS. The addition of a novel MMO loss component, which forces unimodal embeddings to provide independent and complementary information to the fused prediction, further improves prognostic performance. Our DOF model incorporating radiology, histology, and genomic data significantly stratifies glioma patients by OS within outcome-associated subsets, offering additional granularity to routine clinical markers. DOF can be applied to any number of cancer domains, modality combinations, or new clinical endpoints including treatment response.

\printbibliography
\pagebreak
\section{Supplemental Information}
\beginsupplement

\begin{figure}[!hbt]
\centering
\includegraphics[width=.8\textwidth]{radfeaturenet_vgg.png}
\caption{Radiology FeatureNet combining images and features from Gd-T1w and T2w-FLAIR scans. }\label{figS1}
\end{figure}

  \begin{figure}[!hbt]
 \centering
\includegraphics[width=\textwidth]{MICCAI_patient_selection_2021.png}
\caption{Patient selection flowchart.} \label{figS2}
\end{figure}

\begin{table}[!hbt]
\caption{List of handcrafted radiology features. }
\label{tab:tabs1}
\resizebox{\textwidth}{!}{%
\begin{tabular}{@{}lll@{}}
\toprule
\textbf{Feature name/number} & \textbf{Feature Description}                    & \textbf{Summarization of annotated regions} \\ \midrule
No. regions (f1, f2)         & \# annotated lesions on Gd-T1w, edema on T2w-FLAIR      & N/A                                         \\
Volume (f3-f8)               & Volume of 3D ROI, measured in mm$^3$              & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\
Longest axis (f9-f14)        & Longest distance between a contour’s vertices   & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\
SA/V Ratio (f15-f20)         & Ratio of the surface area to volume.            & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\
Sphericity (f21-f26)         & How closely a region’s shape resembles a sphere & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\
Mean I (f27-f32)             & Mean intensity in contoured region              & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\
10th percentile (f33-f38) & 10th \%  of intensities in contoured region & sum, largest, \& avg on Gd-T1w and T2w-FLAIR \\
90th percentile (f39-f44)    & 90th \% of intensities in contoured region      & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\
Skewness (f45-f50)           & Skewness of intensities in contoured region     & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\
Variance (f51-f56)           & Variance of intensities in contoured region     & sum, largest, \& avg on Gd-T1w and T2w-FLAIR            \\ \bottomrule
\end{tabular}%
}
\end{table}

\begin{table}[]
\tiny
\caption{List of molecular features}
\label{tab:tabs2}
\resizebox{\textwidth}{!}{%
\begin{tabularx}{\textwidth}{lXl}
\toprule
\textbf{Category} & \textbf{Variable}        & \textbf{Value Type}                                \\
\midrule
Gene-level CNV (f1-f41) &
  EGFR, MDM4, MYC, BRAF, EZH2, MET, SMO, KIAA1549, CREB3L2, NTRK1, PRCC, BLM, NTRK3, CRTC3, CDKN2A, CDKN2B, FGFR2, TSHR, TCL1A, TRIP11, GOLGA5, GPHN, DICER1, TCL6, EBF1, ITK, RPL22, CDKN2C, LCP1, RB1, IRF4, FGFR1OP, MLLT4, MYB, ROS1, TNFAIP3, GOPC, CARD11, JAK2, STK11, PTEN  & Categorical\\
Arm-level CNV (f42-f78) &
  1q, 2p, 2q, 3p, 3q, 4p, 4q, 5p, 5q, 6p, 6q, 7p, 7q, 8p, 8q, 9p, 9q, 10p, 10q, 11p, 11q, 12p, 12q, 13q, 14q, 15q, 16p, 16q, 17p, 17q, 18p, 18q, 19p, 20p, 20q, 21q, 22q & Continuous \\
Biomarkers (f79, f80)     & IDH Mutation, 1p/19q Codeletion & Binary \\
\bottomrule
\end{tabularx}%
}
\end{table}

\begin{table}[!hbt]
\caption{List of clinical features. }
\label{tab:tabs3}
\resizebox{\textwidth}{!}{%
\begin{tabular}{@{}ll@{}}
\toprule
\textbf{Variable}                                      & \textbf{Value Type} \\ \midrule
Age (f1)                                               & Continuous          \\
Karnofsky Performance Score (f2)                       & Continuous          \\
Grade (f3)                                             & Categorical         \\
Sex: Male vs. Female (f4)                              & Binary              \\
Treatment: any (f5), radiation (f6), chemotherapy (f7) & Binary              \\
Histological diagnosis: LGG (f8), Astrocytoma (f9), Glioblastoma (f10),  Oligoastrocytoma (f11), Oligodendroglioma (f12) & Binary \\
Race/ethnicity: White vs. Non-white (f13), Hispanic vs. Non-hispanic (f14)                                               & Binary \\ \bottomrule
\end{tabular}%
}
\end{table}

\begin{table}[]
\small
\centering
\caption{Median C-index for fusion models at various MMO loss weightings.}
\label{tab:tabs4}
\begin{tabular}{@{}ccc@{}}
\toprule
\textbf{$\gamma$} & \textbf{  Rad + Path + Gen  }    & \textbf{  Rad + Path + Gen + Clin  } \\ \midrule
0                     & 0.764 ± 0.062 & 0.775 ± 0.061   \\
.1                    & 0.768 ± 0.064 & 0.745 ± 0.068   \\
.25                   & 0.777 ± 0.066 & 0.782 ± 0.066   \\
.5                    & 0.788 ± 0.067 & 0.785 ± 0.077   \\
1                     & 0.779 ± 0.070 & 0.776 ± 0.075   \\
2.5                   & 0.781 ± 0.073 & 0.760 ± 0.072   \\ \bottomrule
\end{tabular}
\end{table}

\begin{table}[]
\centering
\small
\caption{Ablation study investigating impact of components of fusion module for best-performing modality combination (Rad + Path + Gen)}
\label{tab:tabs5}
\resizebox{7.25cm}{!}{%
\begin{tabular}{@{}lll@{}}
\toprule
\textbf{  Attention Gating  } & \textbf{  Combination Strategy  } & \textbf{  Median C-index  } \\
\midrule
Yes              & Tensor Fusion        & 0.79 ± 0.07    \\
No               & Tensor Fusion        & 0.77 ± 0.08    \\
Yes              & Concatenation        & 0.78 ± 0.07    \\
No               & Concatenation        & 0.76 ± 0.07   \\

\end{tabular}%
}
\end{table}

\begin{table}[]

\begin{minipage}{\linewidth}
\tiny
\caption{Correlation-based deep fusion framework \cite{cheerla_deep_2019}, adapted to include radiology. }
\label{tab:tabs6}
\resizebox{\textwidth}{!}{%
\begin{tabularx}{\textwidth}{lX}
\toprule
\textbf{Module} & \textbf{Description}                                                                                \\ \midrule
DNA                                                   & Fully connected (FC) layer. Unmodified from \cite{cheerla_deep_2019}.                                                                                \\
Pathology                                            & Squeezenet applied to histology ROIs, followed by a FC layer. While \cite{cheerla_deep_2019} trained from scratch, we found better results when using pretrained ImageNet weights and freezing the convolutional layers                                    \\
Radiology                                             & Not included in \cite{cheerla_deep_2019}. We applied pre-trained, frozen squeezenet to Gd-T1w and T2w-FLAIR ROIs. FC layers were applied to CNN-extracted features and radiomics features, yielding 128 features from each. These were concatenated and processed by another FC layer \\
Fusion                                          & The output of each unimodal arm is a feature vector of length 256, which were averaged together. The averaged feature representation is then processed by a 10-layer highway network. Unmodified from \cite{cheerla_deep_2019}                                                                                \\
Loss                                                     & Combination of Cox proportional hazard-based loss for prognosis prediction and similarity loss that enforces correlated representations between modalities. Unmodified from \cite{cheerla_deep_2019}.                                                                                \\ \bottomrule

\end{tabularx}%
}

\par
\tiny
*All changes from \cite{cheerla_deep_2019} made by us are noted specifically above. Any further differences from the description of \cite{cheerla_deep_2019} are due to discrepancies between the paper and its codebase.\end{minipage}

\end{table}
\end{document}