# Equivalence of optimal $L^{1}$-inequalities on Riemannian Manifolds 1112010
Mathematics Subject Classification: 35A09, 35B44 222Key words: sharp Sobolev
inequalities, best constant, extremal maps

Jurandir Ceccon 333E-mail addresses<EMAIL_ADDRESS>(J. Ceccon)
Departamento de Matemática, Universidade Federal do Paraná,
Caixa Postal 019081, 81531-990, Curitiba, PR, Brazil
Leandro Cioletti 444E-mail addresses<EMAIL_ADDRESS>(L. Cioletti)
Departamento de Matemática, Universidade de Brasília,
UnB, 70910-900, Brasília, Brazil

Abstract

Let $(M,g)$ be a smooth compact Riemannian manifold of dimension $n\geq 2$.
This paper concerns to the validity of the optimal Riemannian $L^{1}$-Entropy
inequality

${\bf Ent}_{dv_{g}}(u)\leq n\log\left(A_{opt}\|Du\|_{BV(M)}+B_{opt}\right)$

for all $u\in BV(M)$ with $\|u\|_{L^{1}(M)}=1$ and existence of extremal
functions. In particular, we prove that this optimal inequality is equivalent
a optimal $L^{1}$-Sobolev inequality obtained by Druet [7].

##  1 Introduction

The isoperimetric problem on the Euclidean space $\mathbb{R}^{n}$ consist in
finding among all the domains with a given fixed volume one that has the
lowest surface area. The solution in this case is a sphere. This property is
precisely expressed in terms of the Isoperimetric inequality for domains in
$\mathbb{R}^{n}$, that is, if $\Omega$ is a domain with volume $|\Omega|$ and
surface area $|\partial\Omega|$ then

$\frac{|\partial\Omega|}{|\Omega|^{\frac{n-1}{n}}}\geq\frac{1}{K(n,1)}\;,$ (1)

where $K(n,1)^{-1}=n^{(1-1/n)}(\omega_{n-1})^{\frac{1}{n}}$ and $\omega_{n-1}$
denotes the volume of the unit ball in the Euclidean space $\mathbb{R}^{n-1}$.
The equality is attained iff $\Omega$ is a sphere. Observe that the last
statement implies that $K(n,1)$ is the best constant for the inequality (1).
The Isoperimetric inequality shows up on other branchs of mathematics. For
example, it has been used to prove the non-uniqueness of the Gibbs measures of
the Ising model on the lattice
$\mathbb{Z}^{n}\equiv\mathbb{Z}\times\ldots\times\mathbb{Z}$. In this context,
the inequality is generalized as follows. We consider $\mathbb{Z}^{n}$ as a
metric space, where the distance between $x,y\in\mathbb{Z}^{n}$ is defined in
terms of the $\ell_{1}$ norm. We also look at $\mathbb{Z}^{n}$ as a graph
$(\mathbb{Z}^{n},\mathbb{E}^{n})$, where $\mathbb{Z}^{n}$ is the vertex set
and
$\mathbb{E}^{n}\equiv\\{\\{x,y\\}\in\mathbb{Z}^{n}\times\mathbb{Z}^{n};|x-y|_{1}=1\\}$
is the edge set. The discrete Isoperimetric inequality says that for any fixed
integer $n\geq 2$ and any finite subset $\Omega\subset\mathbb{Z}^{n}$, we have
that

$\frac{|\partial\Omega|}{|\Omega|^{\frac{n-1}{n}}}\geq\frac{1}{2n},$

where $\partial\Omega=\\{\\{i,j\\}\in\
\mathbb{E}^{d}:i\in\Omega,j\in\Omega^{c}\\}$ and $|\Omega|$ and
$|\partial\Omega|$ denotes the cardinality of $\Omega$ and $\partial\Omega$,
respectively. Note that $2n$ is exactly the volume of the unit sphere on the
$\ell_{1}$ norm. It is worth pointing out that the proof of this discrete
inequality is similar on spirit of our proof on the continuous setting and is
based on the entropy inequalities. It is possible that the equivalences
obtained here can be extended to the discrete case but we will not develop
this point here. More information about this connection can be found in [16].
For more details about the Isoperimetric inequality see the excellent
Osserman’s work [15] and references therein.

Now we consider a more analytic context. We say that a function $u\in
L^{1}(\mathbb{R}^{n})$ has bounded variation if

$\|Du\|_{BV(\mathbb{R}^{n})}=\sup\left\\{\int_{\mathbb{R}^{n}}u\,\mbox{div}(\varphi)\,dx;\,\varphi\in
C_{0}^{1}(\mathbb{R}^{n},\mathbb{R}^{n}),|\varphi|\leq 1\right\\}<\infty\;.$

The space of all bounded variation functions is denoted by
$BV(\mathbb{R}^{n})$. The optimal Euclidean $L^{1}$-Sobolev inequality in
$BV(\mathbb{R}^{n})$ states that for all $u\in BV(\mathbb{R}^{n})$ we have

$\|u\|_{L^{1^{*}}(\mathbb{R}^{n})}\leq K_{0}\|Du\|_{BV(\mathbb{R}^{n})}\;,$
(2)

where $1^{*}=\frac{n}{n-1}$ is the critical Sobolev exponent and

$K_{0}^{-1}=\inf\\{\|Du\|_{BV(\mathbb{R}^{n})};u\in
BV(\mathbb{R}^{n}),\|u\|_{L^{1^{*}}(\mathbb{R}^{n})}=1\\}\;$

is the best constant for this inequality. This inequality it was studied by
Federer and Fleming [9], Fleming and Rishel [10] and Maz′ja [14]. In this
case, the characteristic functions of the balls are extremal functions for the
optimal $L^{1}$-Sobolev inequality and explicit value for best constant is
given by

$K_{0}=K(n,1)\;.$

This inequalities gains in interest if we realize that the geometric
inequality (1) and the analytic inequality (2) are equivalents. This relation
was pointed out independently by Federer and Fleming [9] and Maz′ja [14].
Let’s move on to the manifold setting.

Let $(M,g)$ be a smooth compact Riemannian manifold without boundary of
dimension $n\geq 2$. We say that $u\in L^{1}(M)$ is a bounded variation
function if

$\|D_{g}u\|_{BV}=\sup\left\\{\int_{M}u\,\mbox{div}_{g}(X)\,dv_{g};X\in\Gamma(TM),|X|_{g}(x)\leq
1\hskip 5.69046pt\mbox{for all}\hskip 5.69046ptx\in M\right\\}<\infty\;,$

where $\Gamma(TM)$ is the set of all vector fields over $M$ with divergent in
$n$-fold Cartesian product $C^{1}(M)\times\cdots\times C^{1}(M)$. We denote by
$BV(M)$ the space of bounded variation functions. The optimal Riemannian
$L^{1}$-Sobolev inequality was obtained by Druet [7]. He proved that for all
$u\in BV(M)$ the following inequality holds

$\|u\|_{L^{1^{*}}(M)}\leq K(n,1)\|D_{g}u\|_{BV(M)}+B(1)\|u\|_{L^{1}(M)}\;,$
(3)

where $K(n,1)$ is the first best constant and $B(1)$ is the second best
constant for the optimal Riemannian $L^{1}$-Sobolev inequality.

Moreover, Druet also proved that the only extremal functions for (3) are the
characteristic function of some $\Omega_{0}\in\Sigma$, where
$\Sigma=\\{\Omega\subset M;\chi_{\Omega}\in BV(M)\\}\;.$ Finally, Druet
observed that for $\Omega\in\Sigma$, the optimal inequality (3) is equivalent
to the following isoperimetric inequality

$|\Omega|^{\frac{n-1}{n}}\leq K(n,1)|\partial\Omega|+B(1)|\Omega|,$ (4)

where $|\Omega|$ and $|\partial\Omega|$ are the Riemannian volume and area of
$\Omega$ and $\partial\Omega$, respectively and $B(1)$ is the best constant in
this inequality. For more details about this equivalence and additional
references see Druet [7].

The aim of this work is to show that these inequalities are also equivalent to
the Entropy and Gagliardo-Nirenberg inequalities on both context Euclidean and
Riemannian, and determine the equivalence among its extremal functions.

##  2 Euclidean inequality

For any function $u\in BV(\mathbb{R}^{n})$ we have that its entropy, with
respect to the Lebesgue measure, is well defined and given by the following
expression

${\bf Ent}_{dx}(u)=\int_{\mathbb{R}^{n}}|u|\log|u|dx\;.$

The optimal Euclidean entropy inequality states that for all $u\in
BV(\mathbb{R}^{n})$ with $\|u\|_{L^{1}(\mathbb{R}^{n})}=1$,

${\bf Ent}_{dx}(u)\leq n\log\left(L(n,1)\|\nabla
u\|_{BV(\mathbb{R}^{n})}\right)\;,$ (5)

where $L(n,1)$ is the best constant for this inequality. Note that the
existence of this best constant is guaranteed by the same argument we use in
(2). This inequality it was studied by Beckner [5] and Ledoux [13]. For the
optimal entropy inequality, Beckner showed that the extremal functions for (5)
are the characteristic functions of the Euclidean balls
$B(x_{0},r)=\\{x\in\mathbb{R}^{n};\|x-x_{0}\|<r\\}$, with $|B(x_{0},r)|=1$. So
in this case, it is easy to check that

$K(n,1)=L(n,1)\;.$ (6)

In addition, if $\Omega$ is a domain in $\mathbb{R}^{n}$ and $\lambda>0$ is
chosen so that

$\lambda\int_{M}\chi_{\Omega}dx=1$

then by using $\lambda\chi_{\Omega}$ in (5), we obtain the Isoperimetric
inequality (1).

Now we show that the entropy inequality can be obtained by the Gagliardo-
Nirenberg inequality. Let $1\leq q<r<1^{*}$. By using the interpolation
inequality and the optimal Euclidean $L^{1}$-Sobolev inequality, we obtain the
Euclidean Gagliardo-Nirenberg inequality, which says that for any function
$u\in BV(\mathbb{R}^{n})$ we have

$\|u\|_{L^{r}(\mathbb{R}^{n})}^{\frac{1}{\theta}}\leq
K(n,1)\|Du\|_{BV(\mathbb{R}^{n})}\|u\|_{L^{q}(\mathbb{R}^{n})}^{\frac{1-\theta}{\theta}}\,,$
(7)

where $\theta=\frac{n(r-q)}{r(q(1-n)+n)}\in(0,1)$. Let

$A(n,q,r)^{-1}=\inf\left\\{\|Du\|_{BV(\mathbb{R}^{n})}\|u\|_{L^{q}(\mathbb{R}^{n})}^{\frac{1-\theta}{\theta}};u\in
BV(\mathbb{R}^{n}),\|u\|_{L^{r}(\mathbb{R}^{n})}=1\right\\}$

be the best possible constant for this inequality. If $\chi_{B(0,r)}$ denotes
the characteristic function of the Euclidean ball of radius $r>0$, then we
have that $\chi_{B(0,r)}\in BV(\mathbb{R}^{n})$ and an easy computation shows
that

$\|\chi_{B(0,r)}\|_{L^{r}(\mathbb{R}^{n})}^{\frac{1}{\theta}}=K(n,1)\|D\chi_{B(0,r)}\|_{BV(\mathbb{R}^{n})}\cdot\|\chi_{B(0,r)}\|_{L^{q}(\mathbb{R}^{n})}^{\frac{1-\theta}{\theta}}\;.$

Note that the above equality actually proved that

$A(n,q,r)=K(n,1)\;,$ (8)

for all $1\leq q<r<1^{*}$. Therefore the inequality (7) is in fact the optimal
Euclidean $L^{1}$-Gagliardo-Nirenberg inequality and characteristic function
of the balls are extremal functions for (7).

Proceeding, with minor modifications, as in [4] (see also Section 3) we can
verify that the optimal Euclidean $L^{1}$-Gagliardo-Nirenberg inequality
implies the optimal Euclidean $L^{1}$-entropy inequality (5). Piecing together
these information we conclude that (1) $\Rightarrow$ (2) $\Rightarrow$ (7)
$\Rightarrow$ (5) $\Rightarrow$ (1).

Before proceed, we remark that the equality (8) is the key point in Section 3
to prove the main result of this paper.

##  3 The Riemannian inequality and the main result

Using the interpolation inequality and the optimal Riemannian $L^{1}$-Sobolev
inequality we get for any function $u\in BV(M)$, $1\leq q<r<1^{*}$ that

$\|u\|_{L^{r}(M)}^{\frac{1}{\theta}}\leq\left(K(n,1)\|D_{g}u\|_{BV(M)}+B(1)\|u\|_{L^{1}(M)}\right)\|u\|_{L^{q}(M)}^{\frac{1-\theta}{\theta}}\,,$
($I_{q,r}(K(n,1),B(1))$)

where $\theta=\frac{n(r-q)}{r(q(1-n)+n)}$ is the interpolation parameter. The
first Riemannian $L^{1}$-Gagliardo-Nirenberg best constant is defined by

$A_{opt}=\inf\\{A\in\mathbb{R}:\;\mbox{there exists}\hskip
5.12128ptB\in\mathbb{R}\hskip 5.12128pt\mbox{such that}\hskip
5.12128ptI_{q,r}(A,B)\hskip 5.12128pt\mbox{is valid}\\}\,.$

Using a partition unity and a similar argument as in [8] together with (8), we
can verify that

${\cal A}_{opt}=A(n,q,r)=K(n,1)\;,$

for all $1\leq q<r<1^{*}$. So this equality shows that the first optimal
Riemannian $L^{1}$-Gagliardo-Nirenberg inequality ($I_{q,r}(K(n,1),B(1))$) is
valid for all $u\in BV(M)$ and $K(n,1)$ is the first best constant.

It follows from [7] that every extremal function for the Sobolev inequality is
of the form $u_{0}=\lambda\chi_{\Omega_{0}}\;,$ for some
$\lambda\in\mathbb{R}$ and $\Omega_{0}\in\Sigma$. We see at once that for such
functions the equality in $I_{q,r}(K(n,1),B(1))$ is verified. Consequently the
second Riemannian $L^{1}$-Gagliardo-Nirenberg best constant is given by

$B(1)=\inf\\{B\in\mathbb{R};I_{q,r}(K(n,1),B)\hskip 5.69046pt\mbox{is
valid}\\}\;.$

So we have that optimal Riemannian $L^{1}$-Gagliardo-Nirenberg inequality

$\|u\|_{L^{r}(M)}^{\frac{1}{\theta}}\leq\left(K(n,1)\|D_{g}u\|_{BV(M)}+B(1)\|u\|_{L^{1}(M)}\right)\|u\|_{L^{q}(M)}^{\frac{1-\theta}{\theta}}\,,$
(9)

is valid for all $u\in BV(M)$ and $1\leq q<r<1^{*}$. Notice that the Sobolev
extremal functions are also Gagliardo-Nirenberg extremal functions.

In the sequel we show how to use the inequality (9) (with $q=1$) to obtain the
optimal Riemannian $L^{1}$-entropy inequality. The proof is based on the
Bakry, Coulhon, Ledoux and Sallof-Coste argument given in [4]. Consider the
optimal Riemannian $L^{1}$-Gagliardo-Nirenberg inequality

$\|u\|_{L^{r}(M)}^{\frac{1}{\theta}}\left(\leq
K(n,1)\|D_{g}u\|_{BV(M)}+B(1)\|u\|_{L^{1}(M)}\right)\|u\|_{L^{1}(M)}^{\frac{1-\theta}{\theta}}\;,$

for all $u\in BV(M)$. By taking the logarithm on both sides above and use the
definition of $\theta$, we get

$\frac{r(1-n+n)}{n}\frac{1}{r-1}\log\left(\frac{\|u\|_{L^{r}(M)}}{\|u\|_{L^{1}(M)}}\right)\leq\log\left(K(n,1)\frac{\|D_{g}u\|_{BV(M)}}{\|u\|_{L^{1}(M)}}+B(1)\right).$

Taking the limit when $r\to 1^{+}$, on the above expression, we obtain

$\frac{1}{n}\lim_{r\rightarrow
1^{+}}\frac{1}{r-1}\log\left(\frac{\|u\|_{L^{r}(\mathbb{R}^{n})}}{\|u\|_{L^{1}(\mathbb{R}^{n})}}\right)\leq\log\left(K(n,1)\frac{\|D_{g}u\|_{L(M)}}{\|u\|_{L^{1}(M)}}+B(1)\right)\;.$

To evaluate the remainder limit, we first observe that

$\displaystyle\log\left(\frac{\|u\|_{L^{r}(M)}}{\|u\|_{L^{1}(M)}}\right)$
$\displaystyle=\frac{1}{r}\log(\|u\|_{L^{r}(M)}^{r})-\log(\|u\|_{L^{1}(M)})$
$\displaystyle=\frac{1-r}{r}\log(\|u\|_{L^{1}(M)})+\frac{1}{r}\left(\log(\|u\|_{L^{r}(M)}^{r})-\log(\|u\|_{L^{1}(M)})\right).$

Next, we apply two times the mean value theorem, obtaining

$\lim_{r\rightarrow
1^{+}}\frac{1}{r-1}\log\left(\frac{\|u\|_{L^{r}(M)}}{\|u\|_{L^{1}(M)}}\right)=\int_{M}\frac{|u|}{\|u\|_{L^{1}(M)}}\log\left(\frac{|u|}{\|u\|_{L^{1}(M)}}\right)dx.$

From the above equation it follows that

$\int_{M}|u|\log(|u|)dv_{g}\leq
n\log\left(K(n,1)\|D_{g}u\|_{BV(M)}+B(1)\right)\ ,$

for all $u\in BV(M)$ with $\|u\|_{L^{1}(M)}=1$. As in the previous section we
define

${\bf Ent}_{dv_{g}}(u)=\int_{M}|u|\log|u|dv_{g}\;.$

As a consequence of the previous inequality we have that

${\bf Ent}_{dv_{g}}(u)\leq n\log\left(L(n,1)\|D_{g}u\|_{BV(M)}+B(1)\right)\;,$
($Ent(L(n,1),B(1)$)

for all $u\in BV(M)$ with $\|u\|_{L^{1}(M)}=1$. We shall remember that
$K(n,1)=L(n,1)$. Now we consider the optimal Riemannian $L^{1}$-entropy
inequality

${\bf Ent}_{dv_{g}}(u)\leq n\log\left(L_{opt}\|D_{g}u\|_{BV(M)}+B\right)$

where $u\in BV(M)$ with $\|u\|_{L^{1}(M)}=1$, $B\in\mathbb{R}$ and the first
Riemannian $L^{1}$-entropy best constant is defined by

$L_{opt}=\inf\\{A\in\mathbb{R}:\;\mbox{there exists}\hskip
5.12128ptB\in\mathbb{R}\hskip 5.12128pt\mbox{such that}\hskip
5.12128ptEnt(A,B)\hskip 5.12128pt\mbox{is valid}\\}\,.$

From the definition of best constant and the validity of $Ent(L(n,1),B(1))$
follows that $L_{opt}\leq L(n,1)$. On the other hand, the equality between
these two constants requires a proof.

Assuming that $Ent(L_{opt},B)$ is valid for some $B\in\mathbb{R}$, one can
define the second Riemannian $L^{1}$-Entropy best constant by

$B_{opt}=\inf\\{B\in\mathbb{R};Ent(L_{opt},B)\hskip 5.69046pt\mbox{is
valid}\\}\;.$

Our main result states that $Ent(L_{opt},B_{opt})$ is valid, moreover the
Riemannian first best constant is equals to the Euclidean first best constant.
We also prove that the second best constant for Sobolev and entropy
inequalities are the same.

######  Theorem 1.

Let $(M,g)$ be a smooth compact Riemannian manifold without boundary of
dimension $n\geq 2$. Then $Ent(L_{opt},B(1))$ is valid. In addition,
$L_{opt}=L(n,1)$ and $B_{opt}=B(1)$.

Remark. Together with the results in [7] the Theorem 1 can be used to complete
the equivalence of the four $L^{1}$ optimal inequalities considered on this
work.

Indeed, in [7] we have that isoperimetric $\Rightarrow$ Sobolev. In the
Section 3 we shown that Sobolev $\Rightarrow$ Gagliardo-Nirenberg and
Gagliardo-Nirenberg + Theorem 1 $\Rightarrow$ entropy. Finally, by taking
$\Omega\in\Sigma$ and $\lambda>0$ such that

$\int_{M}\lambda\;\chi_{\Omega}dv_{g}=1\;,$

and replacing $\lambda\;\chi_{\Omega}$ in $Ent(L(n,1),B(1))$ we get the
optimal isoperimetric inequality.

So the four $L^{1}$-optimal inequalities considered here are all equivalent
and we have that (4) $\Rightarrow$ (3) $\Rightarrow$ ($I_{q,r}(K(n,1),B(1))$)
$\Rightarrow$ ($Ent(L(n,1),B(1)$) $\Rightarrow$ (4) .

Proof of Theorem 1: We proceed to show that

$L_{opt}=L(n,1)\;.$

From this equality it may be concluded that $Ent(L(n,1),B(1))$ is the best
$L^{1}$-entropy. To prove the above equality, we use the definition of
$L_{opt}$. That is, by definition, for all $s>0$ there exists $B(s)$ such that

$\int_{M}|u|\log(|u|)dv_{g}\leq
n\log\left((L_{opt}+s)\|D_{g}u\|_{BV(M)}+B(s)\right)\ ,$

for all $u\in BV(M)$ with $\|u\|_{L^{1}(M)}=1$. This is equivalent to the
following inequality

$\frac{1}{\|u\|_{L^{1}(M)}}\int_{M}|u|\log(|u|)dv_{g}+(n-1)\log(\|u\|_{L^{1}(M)})\leq
n\log\left((L_{opt}+s)\|D_{g}u\|_{BV(M)}+B(s)\int_{M}|u|dv_{g}\right)$ (10)

for all $u\in BV(M)$. Let $r>0$ be such that

$\int_{B(0,r)}dx=1\;,$

where $B(0,r)=\\{x\in\mathbb{R}^{n};\|x\|<r\\}$ is the Euclidean ball. For
$\varepsilon>0$, we define
$\chi_{\varepsilon}=(r\varepsilon)^{-n}\chi_{B_{g}(0,\varepsilon)}$, where
$\chi_{B_{g}(0,\varepsilon)}$ is the characteristic function of the Riemannian
ball $B_{g}(0,\varepsilon)=\\{x\in M;d_{g}(x_{0},x)\leq\varepsilon\\}$ and
$x_{0}\in M$. In the normal coordinates in around $x_{0}$, we have

$\int_{M}\chi_{\varepsilon}dv_{g}=(\varepsilon
r)^{-n}\int_{B(x_{0},\varepsilon)}dv_{g}=\int_{B(0,r)}dv_{g_{\varepsilon}}$

and $dv_{g_{\varepsilon}}=\sqrt{(g_{ij}(\exp_{x_{0}}(\varepsilon
x)))}dx\rightarrow dx$ when $\varepsilon\rightarrow 0$, where $g_{ij}$ are the
coefficients of the metric $g$ in the normal coordinates. We also have

$\|D_{g}\chi_{\varepsilon}\|_{BV(M)}=(\varepsilon r)^{-n}\int_{\partial
B(x_{0},\varepsilon)}dv_{g}=(\varepsilon r)^{-1}\int_{\partial
B(0,r)}dv_{g_{\varepsilon}}\;.$

By replacing this identity in (10) we get

$-n\log(\varepsilon
r)+(n-1)\log|B(0,r)|_{g_{\varepsilon}}\leq-n\log(\varepsilon
r)+n\log\left((L_{opt}+s)|\partial B(0,r)|_{g_{\varepsilon}}+\varepsilon
rB(s)|B(0,r)|_{g_{\varepsilon}}\right)\;.$

Because $|B(0,r)|_{\xi}=1$ ($\xi$ denotes the Euclidean metric) and by taking
the limits when $\varepsilon\rightarrow 0$ and $s\rightarrow 0$ in this order,
we obtain

$0\leq n\log(L_{opt}|\partial B(0,r)|_{\xi})\;.$

By the choice of $r$ and remembering again that the characteristic function of
the Euclidean ball $B(0,r)$ is the extremal function in (5), we have now
reached the following identity

$0=n\log(L(n,1)|\partial B(0,r)|_{\xi})\;.$

From where we conclude that $L(n,1)\leq L_{opt}$, which proves the desired
equality. Thus ($Ent(L(n,1),B(1)$) is the optimal Riemannian $L^{1}$-entropy.

Now we compute the second best constant. Consider $\Omega_{0}\in\Sigma$ and
$k>0$, such that

$\int_{M}k\chi_{\Omega_{0}}dv_{g}=1\;,$

such that, $\chi_{\Omega_{0}}$ is the extremal function in the optimal Sobolev
inequality (2). We can immediately check that the optimal inequality
$Ent(L(n,1),B(1))$ becomes an equality when we evaluate in the function
$k\chi_{\Omega_{0}}$. Therefore $B(1)=B_{opt}$.

##  4 Equivalence between extremal functions

As we already observed, the extremal functions for the optimal $L^{1}$-Sobolev
inequality (3) are characteristic functions. We also remarked that these
functions are also extremal for the optimal $L^{1}$-Gagliardo-Nirenberg
inequality ($I_{q,r}(K(n,1),B(1))$). We can also prove, using the limit
process employed in the Section 3, that the extremal functions for the optimal
$L^{1}$-Gagliardo-Nirenberg ($I_{q,r}(K(n,1),B(1))$) are the extremal
functions for the optimal $L^{1}$-Entropy inequality ($Ent(L(n,1),B(1)$).

If we prove that the extremal functions for the optimal $L^{1}$-Entropy
inequality ($Ent(L(n,1),B(1)$) are the extremal functions for the optimal
$L^{1}$-Sobolev inequality (3), we have that the set of the extremal functions
for the four optimal $L^{1}$ inequalities considered here are the same.

Claim. If $u_{0}$ is an extremal function for ($Ent(L(n,1),B(1)$)
$\Longrightarrow$ $u_{0}$ is an extremal function for (2). In fact, from the
Jensen Inequality for any $u_{0}$ such that $\|u_{0}\|_{L^{1}(M)}=1$ we get

$\log\int_{M}|u_{0}|^{1^{*}}dv_{g}=\log\int_{M}|u_{0}|^{1^{*}-1}|u_{0}|dv_{g}\geq\int_{M}\log(|u_{0}|^{\frac{1}{n-1}})|u_{0}|dv_{g}=\frac{1}{n-1}\int_{M}\log(|u_{0}|)|u_{0}|dv_{g}.$

By using ($Ent(L(n,1),B(1)$) it follows that

$\log\int_{M}|u|^{1^{*}}dv_{g}\geq
1^{*}\log\left(K(n,1)\|Du_{0}\|_{BV(M)}+B(1)\right),$

that is,

$\|u_{0}\|_{L^{1^{*}}(M)}\geq K(n,1)\|Du_{0}\|_{BV(M)}+B(1).$

This inequality shows that $u_{0}$ is an extremal function for (3).

Acknowledgments. The first author was partially supported by CAPES through
INCTmat and second author is supported by FEMAT.

## References

  * [1]
  * [2] R. Adams - General logarithmic Sobolev inequalities and Orlicz embedding, J. Funct. Anal. 34, 292-303 (1979).
  * [3] T. Aubin - Problèmes isopérimétriques et espaces de Sobolev, J. Differential Geom. 11, 573-598 (1976).
  * [4] D. Bakry, T. Coulhon, M. Ledoux, L. Sallof-Coste, Sobolev inequalities in disguise, Indiana J. Math., 44 (4), 1033-1074 (1995).
  * [5] W. Beckner - Geometric asymptotics and the logarithmic Sobolev inequality, Forum Math. 11, 105-137 (1999).
  * [6] E. Carlen - Superadditivity of Fisher’s information and logarithmic Sobolev inequalities, J. Funct. Anal. 101, 194-211 (1991).
  * [7] O. Druet - Isoperimetric Inequalities on Compact Manifolds, Geom. Dedicata 90, 217-236 (2002).
  * [8] O. Druet, E. Hebey, M. Vaugon - Optimal Nash’s inequalities on Riemannian manifolds: the influence of geometry, Int. Math. Res. Not. 14, 735-779 (1999).
  * [9] H. Federer, W. H. Fleming - Normal and integral currents, Ann. of Math. 72, 458-520 (1960).
  * [10] W. H. Fleming, R. Rishel - An integral formula for total gradient variation, Arch. Math. 11, 218-222 (1960).
  * [11] L. Gross - Logarithmic Sobolev inequalities, Amer. J. Math. 97, 1061-1083 (1975).
  * [12] E. Hebey, M. Vaugon - Meilleures constantes dans le théorème d’inclusion de Sobolev, Ann. Inst. H. Poincaré. 13, 57-93 (1996).
  * [13] M. Ledoux - Isoperimetry and Gaussian analysis, Lectures on Probability Theory and Statistics (Saint-Flour, 1994), Lecture Notes in Mathematics, v. 1648, Springer, Berlin, 165-294, (1996).
  * [14] V. G. Maz′ja - Classes of domains and imbedding theorems for function spaces, Soviet Math. Dokl. 1, 882-885 (1960).
  * [15] R. Osserman - The isoperimetric inequality, Bull. Amer. Math. Soc. 84, 1182-1238 (1978).
  * [16] L. Saloff-Coste Sobolev inequalities in familiar and unfamiliar settings, Sobolev spaces in mathematics. I, Springer 299–343 (2009).
  * [17] F. Weissler - Logarithmic Sobolev inequalities for the heat-diffusion semigroup, Trans. Amer. Math. Soc. 237, 255-269 (1978).


# clotFoam: An Open-Source Framework to Simulate Blood Clot Formation Under
Arterial Flow

David Montgomery  Federico Municchi  Karin Leiderman  Department of Applied
Mathematics and Statistics, Colorado School of Mines, 1500 Illinois St,
Golden, CO 80401, United States of America  Department of Mechanical
Engineering, Colorado School of Mines, 1500 Illinois St, Golden, CO 80401,
United States of America  Department of Mathematics, University of North
Carolina at Chapel Hill, 216 Lenoir Dr, Chapel Hill, NC 27599, United States
of America  Computational Medicine Program, University of North Carolina at
Chapel Hill, 216 Lenoir Dr, Chapel Hill, NC 27599, United States of America

###### Abstract

Blood clotting involves the coupled processes of platelet aggregation and
coagulation. Simulating clotting under flow in complex geometries is
challenging due to multiple temporal and spatial scales and high computational
cost. _clotFoam_ is an open-source software developed in OpenFOAM that employs
a continuum model of platelet advection, diffusion, and aggregation in a
dynamic fluid environment and a simplified coagulation model with proteins
that advect, diffuse, and react within the fluid and with wall-bound species
through reactive boundary conditions. Our framework provides the foundation on
which one can build more complex models and perform reliable simulations in
almost any computational domain.

###### keywords:

Blood clotting , Platelet aggregation , Coagulation , Hemostasis , Multiscale
modeling , OpenFOAM

††journal: SoftwareX

## Current code version

Current code version |  1.0
---|---
Permanent link to code/repository used for this code version |  https://github.com/d-montgomery/clotFoam
Code Ocean compute capsule |
Legal Code License |  GNU GPL V3
Code versioning system used |  git
Software code languages, tools, and services used |  C++, C, MPI, GNU Make
Compilation requirements, operating environments & dependencies |  OpenFOAM-v9
If available Link to developer documentation/manual |
Support email for questions |<EMAIL_ADDRESS>
Table 1: Code metadata.

##  1 Motivation and significance

Blood clotting is the body’s response to prevent bleeding from an injured
blood vessel. The clotting process involves two main components: platelet
aggregation and coagulation. Platelet aggregation is a primarily physical
process where platelets adhere to the injured vessel wall, and become
activated by receptors that interact with proteins embedded in the wall.
Activated platelets release agonists, such as adenosine diphosphate (ADP),
which can activate and recruits more platelets to the injury where they begin
to form a platelet plug. Coagulation is a biochemical process involving dozens
of enzymatic reactions that occur in the fluid, on activated platelet
surfaces, and on injured portions of the vessel wall. Coagulation features the
interplay of positive and negative feedback loops that work collectively to
encourage thrombus growth in a self-regulating manner. The reactions culminate
in the generation of the enzyme thrombin on activated platelet surfaces, which
is a strong platelet activation agonist, a key player in positive feedback,
and converts fibrinogen into fibrin, which polymerizes and forms a stabilizing
fibrin gel on the platelet plug. Platelets play a critical role in
coagulation, as thrombin generation and inhibition is strongly regulated by
their activated surfaces [1, 2, 3, 4]. Thrombin generation is often used as a
clinical indicator of healthy clotting, as without thrombin, a clot is
typically leaky and unstable. The interested reader can find more information
about the blood clotting process and previously developed mathematical models
in a number of reviews published elsewhere [5, 6, 7, 8, 9].

Spatio-temporal continuum models that employ computational fluid dynamics
(CFD) have gained widespread usage to study blood clotting in devices [10], in
common microfluidic assays [11, 12], in aneurysms [13] and bleeding [14, 15],
and to understand the effects of flow and transport on the clotting process
overall [3, 16, 17, 12]. A significant proportion of spatial-temporal models
of clotting have been implemented using in-house codes [3, 16, 15], commercial
tools with high licensing costs [11], or are not publicly available [10, 18,
13, 12], posing a challenge for researchers with limited resources. There is
growing interest in open source and freely available software tools such as
OpenFOAM [19]. This platform provides a versatile object-oriented toolkit for
developing and constructing CFD software, featuring a diverse range of solvers
and discretization schemes for general grids and parallel computing.

A few in-house codes of clotting models are freely available [20, 21], but are
quite model-specific and not necessarily made for other researchers to easily
build upon. In this work we present _clotFoam_ , a cell-centered finite volume
solver that provides a flexible framework for simulating and easily extending
a reduced model of blood clotting under flow. The software was developed with
OpenFOAM libraries and can be perceived as an extension of the transient fluid
solver _icoFoam_ , which utilizes the PISO algorithm [22] for decoupling
pressure from the fluid velocity. _clotFoam_ incorporates several additions,
including a Darcy term in the fluid equations, five advection-diffusion-
reaction (ADR) equations that describe platelet aggregation and the release of
ADP from platelet stores upon their activation, as well as twelve ADR
equations that represent a reduced model of platelet surface-mediated
coagulation. The coupling of reactions among the platelet and biochemical
species is managed via object-oriented programming and a modified Runga-Kutta
method, employing both field and patch field operations for the computation of
reactions within the fluid and on the surface of the vessel wall. The
following sections of this manuscript will detail the mathematical model
employed by _clotFoam_ , and provide a comprehensive description of the
software. Two examples of the application of _clotFoam_ in 2-D thrombosis and
3-D hemostasis simulations are presented. Finally, we discuss the potential
implications and benefits of the _clotFoam_ solver for the broader modeling
and simulation community.

##  2 Model description

The _clotFoam_ software is based on our previous model that used a continuum
description of fluid, platelets, and platelet aggregation under flow [3]. The
model describes blood as an incompressible Newtonian fluid that is governed by
the Navier-Stokes-Brinkman equations:

$\displaystyle\rho\partialderivative{\overrightarrow{\mathbf{u}}}{t}+\rho(\overrightarrow{\mathbf{u}}\cdot\gradient)\overrightarrow{\mathbf{u}}$
$\displaystyle=-\gradient\mathrm{p}+\mu\nabla^{2}\overrightarrow{\mathbf{u}}-\mu\alpha(\theta^{B})\overrightarrow{\mathbf{u}},$
(1) $\displaystyle\gradient\cdot\overrightarrow{\mathbf{u}}$
$\displaystyle=0,$ (2)

where $\overrightarrow{\mathbf{u}}(\overrightarrow{\mathbf{x}},t)$ is the
fluid velocity, $\mathrm{p}(\overrightarrow{\mathbf{x}},t)$ is pressure,
$\rho$ is the fluid density, and $\mu$ is the dynamic viscosity. The Darcy
term, $-\alpha(\theta^{B})\overrightarrow{\mathbf{u}}$, represents a
frictional resistance to the fluid caused by a growing mass of bound
platelets. The variable $\theta^{B}$ is the ratio of the sum of bound
platelets to the maximum packing limit $P_{\text{max}}$. The permeability of
the mass of bound platelets, $\alpha(\theta^{B})$, decreases as $\theta^{B}$
increases, as it satisfies the Carman-Kozeny relation
$\alpha(\theta^{B})=C_{CK}(0.6\theta^{B})^{2}/(1-0.6\theta^{B})^{3}$, where
$C_{CK}=10^{6}\text{mm}^{-2}$.

Platelets are modeled as number densities (number per volume), eliminating the
need for tracking individual platelets throughout the simulation. The dynamics
of platelet aggregation are described generally by the following hindered ADR
equation:

$\displaystyle\partialderivative{P^{k}}{t}$
$\displaystyle=-\divergence\\{W(\theta^{T})(\overrightarrow{\mathbf{u}}P^{k}-D_{P}\gradient
P^{k})\\}+S_{k}$ (3)

where $P^{k}(\overrightarrow{\mathbf{x}},t)$ is the $k^{th}$ platelet species
and $S_{k}$ is a source/sink term that accounts for the transitions between
different states of the platelet species. The full system of equations in
_clotFoam_ is in A. To ensure that the number of platelets at a location
$\overrightarrow{\mathbf{x}}$ does not exceed a maximum packing limit,
$P_{\text{max}}$, the platelet size is considered using three phenomolgical
functions: the hindered transport function $W(\theta^{T})$, the adhesion
region $k_{\text{adh}}(\overrightarrow{\mathbf{x}})$, and the binding affinity
function $g(\eta)$. A detailed description of these functions can also be
found in A.

The rates of activation of mobile-unactivated platelets by chemical agonists
ADP and thrombin ($E_{2}$) are assumed to satisfy Hill functions of the form
$A(c)=k_{c}^{\text{pla}}\frac{c}{c^{\ast}+c}$. ADP is secreted by newly
activated platelets over a period of 1-5 seconds after activation. The molar
concentration of ADP satisfies:

$\partialderivative{[\text{ADP}]}{t}=-\divergence\Big\\{\overrightarrow{\mathbf{u}}[\text{ADP}]-D_{\text{ADP}}\gradient[\text{ADP}]\Big{\\}}+\sigma_{\text{release}},$
(4)

where $D_{\text{ADP}}$ is the diffusion coefficient, and the source term is
defined as:

$\displaystyle\sigma_{\text{release}}(\overrightarrow{\mathbf{x}},t)$
$\displaystyle=\int_{0}^{\infty}\hat{A}R(\tau)\partialderivative{t}(P^{b,a}+P^{se,a})(\overrightarrow{\mathbf{x}},t-\tau)d\tau.$
(5)

$\hat{A}$ is the total concentration of ADP released by an activated platelet,
and $\hat{A}R(\tau)$ is the rate of release of ADP $\tau$ seconds after
activation. The rate function, $R(\tau)$, utilized by _clotFoam_ is similar to
the one in our previous work [3], but here uses a bell curve centered at 3
seconds, $R(\tau)=\frac{1}{\sqrt{\pi}}\exp(-(\tau-3)^{2})$, and is normalized
such that $\int_{0}^{\infty}R(\tau)\,d\tau=1.$

The reduced coagulation model consists of various biochemical species in molar
concentrations and categorized as fluid-phase, platelet-bound, or
subendothelium-bound. In an effort to enhance the adaptability of _clotFoam_ ,
the number of biochemical species has been reduced from 50 [3, 16] to 12. This
model is an extension of a previously published ODE model [1] that includes
positive feedback and enzyme inhibition, two features necessary to capture the
bursting thrombin behavior observed in coagulation. We extended the ODE model
to a PDE model where reactions occur on two surfaces (subendothelium and
activated platelets) instead of one, and the species are subjected to flow.
The model is detailed in our B, and a schematic is presented in Figure 1. The
reduced model is summarized as follows where $S$ denotes a substrate or
zymogen, and $E$ denotes an enzyme:

  1. 1.

Fluid-phase substrate, $S_{1}$, comes into contact with enzyme, $E_{0}$, bound
to the subendothelium. $E_{0}$ converts the substrate to an enzyme, $E_{1}$.

  2. 2.

$E_{1}$ binds to the activated platelet surface and becomes the bound species,
$E_{1}^{b}$.

  3. 3.

Additional fluid-phase substrates, $S_{2}$, bind to the platelet surface and
become $S_{2}^{b}$. Upon activation by $E_{1}^{b}$, the platelet-bound
substrates, $S_{2}^{b}$, are converted into a second platelet-bound enzyme,
$E_{2}^{b}$, which we consider to be similar to thrombin.

  4. 4.

$E_{2}^{b}$ activates platelet-bound substrates, forming more enzymes in a
positive feedback loop.

  5. 5.

$E_{2}$ activates mobile unactivated platelets.

Figure 1: Schematic representation of the reduced model of thrombin generation
with positive feedback. The light blue circles represent a substrate, yellow
circles denote enzymes, and the red circles represent binding sites on the
platelet surface. Dotted lines indicate transport by the fluid, while the
solid lines represent binding interactions. Note: Platelets have diameters
near 3 $\mu$m, while endothelial cells are typically 50-70 $\mu$m in length
and 0.1-10 $\mu$m in thickness [23].

##  3 Software description

_clotFoam_ is an open-source software distributed under the GNU General Public
License, compiled using the OpenFOAM-v9 libraries. The software can be
compiled on any system where the OpenFOAM-v9 libraries are installed, and full
installation instructions are provided in the repository. The code is written
in C++ and can be adapted to simulate clotting in a wide variety of domains
with few limitations to the mesh. The use of object-oriented programming
enables the management of platelets and biochemical species as objects, which
simplifies the implementation of more complex coagulation and clotting models.
The repository contains two illustrative examples, with the expectation that
its range will expand through community contributions and author updates as
the framework develops further. _clotFoam_ is fully parallelizable for high
performance computing (HPC) using the message passing interface (MPI)
framework. The solution algorithm employed by _clotFoam_ is illustrated in
Figure 2 and is described in more detail in the following subsections.

Figure 2: Flow chart of platelet-mediated coagulation solver within OpenFOAM.

###  3.1 Mesh requirements

To properly define the injury region, the mesh needs to satisfy two
requirements. First, the reactive boundary conditions and subendothelium-bound
species are defined exclusively on a patch known as “injuryWalls”. This patch
is constructed to be the wall of an injury block within the domain using the
blockMesh tool, as depicted in Figure 3. Alternatively, users can use the
topoSet tool to define the injuryWalls patch. Second, the mesh at the injury
site must feature cell widths, heights, and depths no greater than a platelet
diameter $P_{\text{diam}}$, which has a default value of 3 $\mu$m. This
requirement is a consequence of how the adhesion region
$H_{\text{adh}}(\overrightarrow{\mathbf{x}})$ is defined.

Figure 3: Example of thrombus growth in a corresponding domain defined with
blockMesh using three blocks and 16 vertices. The injuryWalls patch is defined
by vertices {1, 9, 10, 2}, which is used by _clotFoam_ to determine the
location of the reactive boundary conditions in the coagulation model.

###  3.2 Numerical methods

The fluids solver in _clotFoam_ is built upon OpenFOAM’s transient fluid
solver _icoFoam_ , which implements the PISO algorithm as a predictor and
corrector method. In each time step, the momentum equation is solved once,
followed by multiple pressure and velocity corrections. During the
discretization of the Navier-Stokes-Brinkman equations (1), the Darcy term is
treated implicitly as a source term. This treatment involves utilizing the
bound platelet fraction, $\theta^{B}$, from the previous time step. By
incorporating the bound platelet fraction in this manner, we not only enhance
the stability of the fluids solver but also ensure that the pressure
corrections in the PISO algorithm are influenced by the presence of the porous
media. Consequently, we are able to accurately capture the influence of the
porous media on the fluid flow, resulting in simulations that exhibit improved
reliability and robustness. Further details regarding other discretization
schemes employed in this work can be found in D.

The model accounts for the transport of both platelet and biochemical species
with equations that incorporate advection, diffusion, and reactions (ADR) with
other species. The reactions occur on a smaller time scale than the transport
processes, specifically, the software’s default reaction time-step is half of
the transport time-step. As such, the software employs a fractional-step
method to decouple the transport from the reaction terms in the equations.
Thus, the reaction equations are solved multiple times during each fluid time-
step. The general ADR equation for each species is:

$\partialderivative{C_{i}}{t}=-\divergence(\overrightarrow{\mathbf{v}}C_{i}-D_{c}\gradient
C_{i})+R_{i}(C_{1},C_{2},...,C_{n}),$ (6)

where $C_{i}$ is the $i^{th}$ species with $i=1,\dots,n$,
$\overrightarrow{\mathbf{v}}$ is a fluid velocity (not necessarily
$\overrightarrow{\mathbf{u}}$ from (1)), $D_{c}$ is the diffusion coefficient,
and $R_{i}$ is a reaction term that can depend on multiple species.

A two-step fractional-step method is used to march equation (6) forward in
time for each time step:

  1. 1.

Solve the transport equation with a temporal step size $\Delta t$:

$\partialderivative{C_{i}}{t}=-\divergence(\overrightarrow{\mathbf{v}}C_{i}-D_{c}\gradient
C_{i}).$ (7)

  2. 2.

Update the solution by solving the coupled reaction equations M_rxn times with
a temporal step size $h=\Delta t/M_{\text{rxn}}$:

$\partialderivative{C_{i}}{t}=R_{i}(C_{1},C_{2},\dots,C_{n}).$ (8)

The transport equations are discretized using the finite volume method (FVM)
as discussed in Appendix D,  while the reaction equations are solved using a
coupled fourth-order Runge-Kutta (RK4) method. The parameter M_rxn determines
how many times the reaction equations are solved per fluid time-step, and is
specified in the $FOAM_CASE/constant/inputParameters file. The default value
is $\texttt{M\\_rxn}=2$, however, it should be noted that this parameter is
dependent on the specific problem and may need to be adjusted for flows with
higher wall-shear rates.

In the mobile platelet equations described in A, the flux vector
$\overrightarrow{\mathbf{j}}=\overrightarrow{\mathbf{u}}P-D_{P}\gradient P$ is
scaled by a hindered transport function $W(\theta^{T})$ to limit the transport
of platelets near the growing thrombus. Prior to FVM discretization of the
advective and diffusive fluxes, the total platelet fraction $\theta^{T}$ must
be interpolated to the cell faces. The choice of interpolation method is
determined by the mechanism of transport. By implementing a combination of
interpolation schemes, the flux of platelets into a spatial location is
effectively constrained, ensuring that the maximum value of the sum of all
platelet species at a spatial location remains below or equal to a maximum
packing density $P_{\text{max}}$. To interpolate the total platelet fraction
for the advective flux, $W(\theta^{T})\overrightarrow{\mathbf{u}}P$, a
downwind scheme is used with respect to the fluid velocity
$\overrightarrow{\mathbf{u}}$. This is because the fluid velocity is assumed
to only be hindered by a thrombus that is downstream. Conversely, for the
diffusive flux, $W(\theta^{T})D_{P}\gradient P$, the total platelet fraction
is interpolated using a localMax scheme, as the diffusion rate within the
thrombus is expected to be smaller than the rate outside of the thrombus.

The ADP equation (4) is not solved using the fractional-step method, because
the source term is updated infrequently. When platelets become bound, they
release ADP into the fluid for up to 6 seconds. The secretion of ADP is
modeled by the source term, $\sigma_{\text{release}}$, as defined in equation
(5), and can be restricted to the interval $\tau\in[0,\,6]$ due to the bell-
shaped distribution of $R(\tau)$. However, the computation of
$\sigma_{\text{release}}$ is memory intensive as the number of newly bound
platelets
$\partialderivative{t}(P^{b,a}+P^{se,a})(\overrightarrow{\mathbf{x}},t_{n}-\tau)$
must be stored in memory for up to 6 seconds. To reduce computational cost, a
coarse discretization of $\tau$ is employed to calculate
$\sigma_{\text{release}}$. The number of newly bound platelets are computed
and stored at a user-specified interval $\Delta\tau$. The discretization of
$\sigma_{\text{release}}$ at time $t=t_{n}$ is then implemented using the
trapezoid rule:

$\displaystyle\sigma_{\text{release}}(\overrightarrow{\mathbf{x}},t_{n})$
$\displaystyle=\int_{0}^{\tau_{f}}\hat{A}R(t^{\prime})\partialderivative{t}(P^{b,a}+P^{se,a})(\overrightarrow{\mathbf{x}},t_{n}-t^{\prime})\,dt^{\prime},$
$\displaystyle=\int_{t_{n}-\tau_{f}}^{t_{n}}\hat{A}R(t_{n}-\tau)\partialderivative{t}(P^{b,a}+P^{se,a})(\overrightarrow{\mathbf{x}},\tau)\,d\tau,\hskip
14.22636pt(\text{by substituting $\tau=t_{n}-t^{\prime}$})$
$\displaystyle\approx\sum_{k=0}^{N_{\tau}-1}\hat{A}\,\frac{\Delta\tau}{2}\bigg{\\{}R(t_{n}-\tau_{k})\partialderivative{t}(P^{b,a}+P^{se,a})(\overrightarrow{\mathbf{x}},\tau_{k})+R(t_{n}-\tau_{k+1})\partialderivative{t}(P^{b,a}+P^{se,a})(\overrightarrow{\mathbf{x}},\tau_{k+1})\bigg{\\}},$
(9)

where $\tau_{f}$ and $\Delta\tau$ are defined as sigma_Tf and sigma_dt
respectively in the inputParameters dictionary. Lastly, the number of newly
bound platelets is approximated as:

$\displaystyle\partialderivative{t}(P^{b,a}+P^{se,a})(\overrightarrow{\mathbf{x}},t_{n})$
$\displaystyle\approx
P_{\text{max}}\frac{\theta^{B}_{n}-\theta^{B}_{n-1}}{t_{n}-t_{n-1}}.$ (10)

###  3.3 Managing platelet and biochemical species with polymorphism

Models of blood clotting typically involve multiple platelet and biochemical
species. For instance, the Leiderman-Fogelson model [3] consists of four
platelet species and 50 biochemical species. To address this complexity,
_clotFoam_ has been developed to accommodate models with any number of
platelet and biochemical species. These models are implemented following a
polymorphic approach, with an abstract base class called Species from which
four classes are derived to align with the species defined in the mathematical
model: Species_platelet, Species_seBound, Species_fluidPhase,
Species_pltBound. The inheritance relationship of these classes is depicted in
Figure 4.

Figure 4: Graph of the Species class and derived classes Species_platelet,
Species_seBound, Species_fluidPhase, Species_pltBound. All derived classes
inherit the public member variables and public member functions defined in the
Species class.

The Species object consists of four pointer lists to facilitate the management
of the subspecies fields, including field values (solutions), previous field
values (solutions from intermediate time steps used in the fractional-step
method and RK4 solver), $k$ values (computed for the RK4 solver), and the
argument of the $k$ values (inputs for the $k$ values in the RK4 solver). The
public member functions of the derived classes enable the setting of pointers
to other Species objects and define the reaction functions specific to their
corresponding subspecies. The function updateKs computes and retains the
reaction term, $R_{i}$, for each subspecies, utilizing input parameters passed
from the RK4 method. The initialization of every derived Species object occurs
within the createFields.H file, and the pointers are set in
setSpeciesPointers.H.

###  3.4 Adapting the framework for different clotting models

_clotFoam_ was developed with the aim of promoting scientific research and
reproducibility in the field of hemostasis, thrombosis, and clotting
disorders. Thus, we illustrate the procedure for introducing a new species
into the software. Consider the addition of a fluid-phase species $I$, that
inhibits thrombin indefinitely. This reaction can be written as:

$I+E_{2}\xrightarrow{k_{I}}E_{2,\text{inh}},$ (11)

where $k_{I}$ is the association rate of inhibitor and enzyme, and
$E_{2,\text{inh}}$ is the resulting inhibited thrombin. Although the addition
of two new species, $I$ and $E_{2,\text{inh}}$, is necessary for incorporating
the desired reaction into the software, we will focus our discussion on the
steps involved in including just the inhibitor species as it applies to any
additional species. The new reaction term that is added to the ADR equation
for $I$ is:

$\text{reaction rate:}-k_{I}IE_{2},$ (12)

which fits the form of equation (6), with a reaction term defined by the
right-hand side. The following sequential steps outline the procedure for
adding $I$ to the software:

  1. 1.

Within the case directory $FOAM_CASE:

    1. (a)

In the 0 directory:

      1. i.

Create a new field, such as fluidPhase_I.

    2. (b)

In the constant/inputParameters file:

      1. i.

Update the number of fluid-phase species, num_fluidPhase.

      2. ii.

Incorporate the parameter kI using appropriate units.

  2. 2.

Within the _clotFoam_ software:

    1. (a)

In the chemConstants.H file:

      1. i.

Read in the parameter kI defined in the inputParameters dictionary.

    2. (b)

In the Species_fluidPhase.H class file:

      1. i.

Update the reaction function updateKs to include

$\texttt{k[4] = -kI*I*E2 },$

where I is a reference to the input necessary for the RK4 solver,

const volScalarField& I = input[4],

defined at the beginning of the updateKs function.

      2. ii.

If necessary, update the pointers in the constructor and in the setPtrs
function (not required for this example).

    3. (c)

In the createFields.H file:

      1. i.

Add the fluidPhase_I field to the 4th index of the PtrList called
fluidPhase.field.

    4. (d)

In the setSpeciesPointers.H file:

      1. i.

If required, update the fluidPhase.setPtrs argument corresponding to step
2(b)ii.

##  4 Illustrative Examples

Here, we present two illustrative examples for which the code and instructions
are included in the repository. These examples showcase the capabilities of
the _clotFoam_ solver in simulating blood clotting phenomena. The first
example is a 2D thrombosis case based on previous results published by our
group [3, 16]. The second example is a 3D simulation of hemostasis that
replicates the microfluidic device described in our previous study by Schoeman
et al. [14]. Additional examples of convergence and validation can be found in
the supplementary material.

###  4.1 Thrombosis in a rectangular channel

To verify the reliability of _clotFoam_ in modeling clotting phenomena, we
compare its results with those in our previous work [3, 16] and simulate
clotting in a 240 $\mu$m long by 60 $\mu$m high rectangular channel with an
approximately 90 $\mu$m long adhesive and reactive patch centered on the
bottom wall. The domain is discretized using a uniform mesh of (128$\times$32)
cells, which is divided into three blocks as depicted in Figure 3 . It should
be noted that the coagulation reactions used in _clotFoam_ are a simplified
version of those used in our previous work, and therefore, the outcomes are
not expected to be identical. Nevertheless, we demonstrate that _clotFoam_
produces similar concentrations of bound thrombin, clot sizes and densities,
all on the same timescale as in our previous studies [3, 16]. This is
illustrated in Figure 5, where we show spatial distributions of substrates,
enzymes, and clots formed after 200, 400, and 600 seconds of clotting
activity.

Figure 5: Close-up of 140 $\mu$m long by 35 $\mu$m high region around the
thrombus to view various species (rows) at times 200, 400 and 600 seconds in
the left, middle and right columns, respectively. The first 3 rows show
spatial concentrations (in nM) of $E_{1}$ and $S_{2}$ and the total thrombin
concentration ($E_{2}+E_{2}^{b}$), which was calculated in ParaView during
post-processing. The bottom row shows the resulting thrombus growth via the
bound platelet fraction in a dynamic fluid environment where the fluid
velocity field is depicted by the white arrows.

###  4.2 Hemostasis in a microfluidic device

To demonstrate the versatility of _clotFoam_ in simulating clotting in various
two and three dimentional domains, we have included a case that replicates the
H-shaped microfluidic device used in Schoeman et al. [14] to model hemostasis.
In their experiments, whole blood is introduced into the right “blood
channel,” while a buffer fluid is introduced into the left vertical channel.
The pressure difference between the channels causes blood to flow through the
horizontal “injury channel,” which is coated with tissue factor and collagen
proteins that initialize coagulation and platelet adhesion/aggregation
respectively. Clots build up in the injury channel without restricting blood
flow in the blood channel. Figure 6 shows snapshots of the 3D domain at four
different times during the clotting process. Due to difficulties visualing the
entire clot in the injury channel, we have displayed slices through the injury
channel to highlight the spatial growth of the clot at various points down the
injury channel. The rectangular plots to the left of each 3D domain are an
enlarged view of the slice at the location indicated by the black arrow and
show the clot distribution and its dynamics over time. The thrombus growth
patterns observed in the experimental results presented by Schoeman et al.
[14] show thrombus formation primarily at the front (right) of the injury
channel while the thrombus generated by _clotFoam_ , shown in Figure 6, shows
the buildup of thrombus closer to the back (left) of the injury channel. We
hypothesize that this discrepancy may be attributed to the absence of shear
dependence in the current platelet aggregation model, and this aspect will be
investigated as part of our future work.

Figure 6: Snapshots of a _clotFoam_ simulation of hemostasis in an H-shaped
microfluidic device for time 400, 800, 1200, 1600 seconds. The simulation
shows mobile platelets entering the right vertical channel and being
transported by the fluid (indicated by arrows) through the horizontal injury
channel. As they pass through the injury, they start to aggregate and form a
platelet plug, which is visualized via the bound platelet fraction
$\theta^{B}$ in the five slices within the injury channel.

##  5 Impact and conclusions

Computational models that simulate clotting phenomena have provided a
significant step toward a better understanding of hemostasis, thrombosis, and
clotting disorders. In this work, we have presented _clotFoam_ , an open-
source software for simulating clotting using the computational fluid dynamics
framework OpenFOAM. To demonstrate the reliability of _clotFoam_ , we compared
its outcomes with our previous computational work studying thrombosis in
rectangular channels [3, 16, 14, 15]. Our results show that _clotFoam_
produces a similar clot structure and thrombin concentrations, thereby
verifying its implementation of a reduced coagulation and platelet aggregation
model. Furthermore, we demonstrated the versatility of _clotFoam_ by
simulating clotting in an H-shaped microfluidic device, which was used in
experiments for modeling hemostasis. While the simulations did not replicate
the experiments exactly, they illustrated the potential of _clotFoam_ for
investigating clotting in various microfluidic geometries.

In conclusion, _clotFoam_ offers a reliable and flexible platform for
simulating clotting phenomena. With the ability to manipulate the platelet and
coagulation cascade models, researchers can use _clotFoam_ to investigate
various aspects of thrombus formation and design microfluidic devices for
studying hemostasis. The open-source nature of _clotFoam_ also allows for
community-driven development and improvement of the software, making it an
accessible tool for researchers in the field.

## Conflict of Interest

No conflict of interest exists: We wish to confirm that there are no known
conflicts of interest associated with this publication and there has been no
significant financial support for this work that could have influenced its
outcome.

## Acknowledgements

We would like to thank our collaborators for their support and guidance: Aaron
Fogelson, Suzanne Sindi, Keith Neeves, Dougald Monroe, and Michael Stobb. This
work was, in part, supported by the National Institutes of Health (R01
HL151984), and the National Science Foundation CAREER program (DMS-1848221).

##  Appendix A Platelet Aggregation Model

From our previous work, the number densities of the four platelet species are
denoted by $P^{m,u},\;P^{m,a},\;P^{b,a},\;P^{b,se}$, where the superscripts
represent mobile unactivated, mobile activated, platelet-bound activated, and
subendothelium-bound activated respectively. Platelet aggregation is described
by the following system of ADR equations:

$\displaystyle\partialderivative{P^{m,u}}{t}$
$\displaystyle=-\underbrace{\divergence\\{W(\theta^{T})(\overrightarrow{\mathbf{u}}P^{m,u}-D_{P}\gradient
P^{m,u})\\}}_{\text{Transport by advection and
\leavevmode\ltxml@oqmark@open\textquotedblleft\penalty
10000\hskip-0.0002pt\hskip
0.0002ptdiffusion\textquotedblright\ltxml@oqmark@close{}}}$ (13)
$\displaystyle\;\;\;\;-\underbrace{k_{\text{adh}}H_{\text{adh}}(\overrightarrow{\mathbf{x}})(P_{\text{max}}-P^{b,se})P^{m,u}}_{\text{Adhesion
to
subendothelium}}-\underbrace{\\{A_{1}([\text{ADP}])+A_{2}(E_{2})\\}P^{m,u}}_{\text{Activation
by ADP or thrombin}},$ $\displaystyle\partialderivative{P^{m,a}}{t}$
$\displaystyle=-\divergence\\{W(\theta^{T})(\overrightarrow{\mathbf{u}}P^{m,a}-D_{P}\gradient
P^{m,a})\\}-k_{\text{adh}}(\overrightarrow{\mathbf{x}})(P_{\text{max}}-P^{b,se})P^{m,a}$
(14)
$\displaystyle\;\;\;\;+\\{A_{1}([\text{ADP}])+A_{2}(E_{2})\\}P^{m,u}-\underbrace{k_{\text{coh}}g(\eta)P_{\text{max}}P^{m,a}}_{\text{Cohesion
to bound platelets}},$ $\displaystyle\partialderivative{P^{b,a}}{t}$
$\displaystyle=-k_{\text{adh}}H_{\text{adh}}(\overrightarrow{\mathbf{x}})(P_{\text{max}}-P^{b,se})P^{b,a}+k_{\text{coh}}g(\eta)P_{\text{max}}P^{m,a},$
(15) $\displaystyle\partialderivative{P^{b,se}}{t}$
$\displaystyle=k_{\text{adh}}H_{\text{adh}}(\overrightarrow{\mathbf{x}})(P_{\text{max}}-P^{b,se})(P^{m,u}+P^{m,a}+P^{b,a}).$
(16)

The hindered transport function is a monotonically decreasing function defined
as $W(\theta^{T})=\tanh(\pi(1-\theta^{T}))$, where $\theta^{T}$ is the ratio
of the sum of all platelet species to $P_{\text{max}}$. The function assumes
that the transport of mobile platelets is only slightly hindered until
$\theta^{T}$ approaches approximately 0.5, after which the transport of mobile
platelets is drastically reduced. The adhesion function
$H_{\text{adh}}(\overrightarrow{\mathbf{x}})$ is chosen to be unity or zero,
and defines the region where mobile platelets can stick to specified walls on
the domain boundary.

Platelet-platelet cohesion is modeled through the parameters $k_{\text{coh}}$
and the binding-affinity function $g(\eta)$, where the latter depends on a
non-dimensional virtual substance $\eta$ that is produced by bound platelets.
In _clotFoam_ , the virtual substance is modeled by diffusing the bound
platelet fraction a distance of $L_{\eta}$ over an interval of $\Delta t$
through the use of a discrete diffusion equation:

$\displaystyle\frac{\eta-\eta_{0}}{\Delta t}=D_{\eta}\nabla^{2}\eta,$ (17)

where $\eta_{0}$ is the current known bound platelet fraction,
$\eta_{0}=\theta^{B}_{n}$. The diffusion constant is defined as
$D_{\eta}=\frac{L_{\eta}^{2}}{4\Delta t}$, and therefore (17) is implemented
as:

$\displaystyle\eta-\frac{{L_{\eta}}^{2}}{4}\nabla^{2}{\eta}=\theta^{B}_{n}.$
(18)

The definition of $\eta$ presented above departs somewhat from the
corresponding quantity in our previous work, and is intended to rectify the
overly diffusive approach used in our previous work [3, 16, 15]. The binding
affinity function is defined as:

$g(\eta)=\frac{g_{0}\,(\eta-\eta_{t})^{3}}{\eta_{\ast}^{3}+(\eta-\eta_{t})^{3}},$
(19)

where $\eta_{t}$ is a threshold value for which there is no binding,
$\eta_{\ast}+\eta_{t}$ indicates the value of $\eta$ for which $g(\eta)$
changes rapidly, and
$g_{0}=\frac{\eta_{\ast}^{3}+(1-\eta_{t})^{3}}{(1-\eta_{t})^{3}}$ so that
$g(1)=1$.

##  Appendix B Coagulation model

The reduced coagulation network is modeled by the following set of enzymatic
reactions:

$\displaystyle S_{1}+E_{0}\xrightleftharpoons[k_{C_{0}}^{-}]{k_{C_{0}}^{+}}$
$\displaystyle\;C_{0}\xrightarrow{k_{C_{0}}^{\text{cat}}}E_{0}+E_{1}$ (20)
$\displaystyle S_{1}+P_{1}$
$\displaystyle\xrightleftharpoons[k_{S_{1}}^{-}]{k_{S_{1}}^{+}}S_{1}^{b}$ (21)
$\displaystyle E_{1}+P_{1}$
$\displaystyle\xrightleftharpoons[k_{E_{1}}^{-}]{k_{E_{1}}^{+}}E_{1}^{b}$ (22)
$\displaystyle S_{2}+P_{2}$
$\displaystyle\xrightleftharpoons[k_{S_{2}}^{-}]{k_{S_{2}}^{+}}S_{2}^{b}$ (23)
$\displaystyle E_{2}+P_{2}$
$\displaystyle\xrightleftharpoons[k_{E_{2}}^{-}]{k_{E_{2}}^{+}}E_{2}^{b}$ (24)
$\displaystyle
S_{2}^{b}+E_{1}^{b}\xrightleftharpoons[k_{C_{1}}^{-}]{k_{C_{1}}^{+}}$
$\displaystyle\;C_{1}\xrightarrow{k_{C_{1}}^{\text{cat}}}E_{1}^{b}+E_{2}^{b}$
(25) $\displaystyle
S_{1}^{b}+E_{2}^{b}\xrightleftharpoons[k_{C_{2}}^{-}]{k_{C_{2}}^{+}}$
$\displaystyle\;C_{2}\xrightarrow{k_{C_{2}}^{\text{cat}}}E_{1}^{b}+E_{2}^{b}$
(26)

Using the law of mass action, a system of twelve differential equations is
derived to govern the reactions of the twelve biochemical species present in
the coagulation cascade. These equations are categorized as fluid-phase,
platelet-bound, and subendothelium-bound for organizational purposes. All
concentrations of the biochemical species are measured in units of nmol/mm3,
and platelet densities are converted to concentrations using Avogadro’s
constant $N_{A}=6.02214076\times 10^{14}$ nmol-1, where $N_{A}$ is the number
of molecules in one mole of a substance.

### Fluid-phase:

$\displaystyle\partialderivative{S_{1}}{t}$
$\displaystyle=-\divergence(\overrightarrow{\mathbf{u}}S_{1}-D\gradient
S_{1})$ (27)
$\displaystyle\;\;\;\;-k_{S_{1}}^{+}\Big{\\{}\tfrac{N_{1}}{N_{A}}(P^{b,a}+P^{se,a})-(S_{1}^{b}+E_{1}^{b}+C_{1}+C_{2})\Big{\\}}S_{1}+k_{S_{1}}^{-}S_{1}^{b},$
$\displaystyle\hskip 14.22636pt\text{ with }-D\,\hat{\mathbf{n}}\cdot\gradient
S_{1}\big{|}_{\partial\Omega_{\text{inj}}}=-k_{C_{0}}^{+}S_{1}E_{0}+k_{C_{0}}^{-}C_{0},$
$\displaystyle\partialderivative{S_{2}}{t}$
$\displaystyle=-\divergence(\overrightarrow{\mathbf{u}}S_{2}-D\gradient
S_{2})$ (28)
$\displaystyle\;\;\;\;-k_{S_{2}}^{+}\Big{\\{}\tfrac{N_{2}}{N_{A}}(P^{b,a}+P^{se,a})-(S_{2}^{b}+E_{2}^{b}+C_{1}+C_{2})\Big{\\}}S_{2}+k_{S_{2}}^{-}S_{2}^{b},$
$\displaystyle\partialderivative{E_{1}}{t}$
$\displaystyle=-\divergence(\overrightarrow{\mathbf{u}}E_{1}-D\gradient
E_{1})$ (29)
$\displaystyle\;\;\;\;-k_{E_{1}}^{+}\Big{\\{}\tfrac{N_{1}}{N_{A}}(P^{b,a}+P^{se,a})-(S_{1}^{b}+E_{1}^{b}+C_{1}+C_{2})\Big{\\}}E_{1}+k_{E_{1}}^{-}E_{1}^{b},$
$\displaystyle\hskip 14.22636pt\text{ with }-D\,\hat{\mathbf{n}}\cdot\gradient
E_{1}\big{|}_{\partial\Omega_{\text{inj}}}=k_{C_{0}}^{\text{cat}}C_{0},$
$\displaystyle\partialderivative{E_{2}}{t}$
$\displaystyle=-\divergence(\overrightarrow{\mathbf{u}}E_{2}-D\gradient
E_{2})$ (30)
$\displaystyle\;\;\;\;-k_{E_{2}}^{+}\Big{\\{}\tfrac{N_{2}}{N_{A}}(P^{b,a}+P^{se,a})-(S_{2}^{b}+E_{2}^{b}+C_{1}+C_{2})\Big{\\}}E_{2}+k_{E_{2}}^{-}E_{2}^{b}.$

### Platelet-bound:

$\displaystyle\partialderivative{S_{1}^{b}}{t}$
$\displaystyle=\underbrace{k_{S_{1}}^{+}\Big{\\{}\tfrac{N_{1}}{N_{A}}(P^{b,a}+P^{se,a})-(S_{1}^{b}+E_{1}^{b}+C_{1}+C_{2})\Big{\\}}S_{1}-k_{S_{1}}^{-}S_{1}^{b}}_{\text{Binding
and unbinding on platelet surface}}$ (31)
$\displaystyle\;\;\;\;-k_{C_{2}}^{+}S_{1}^{b}E_{2}^{b}+k_{C_{2}}^{-}C_{2},$
$\displaystyle\partialderivative{S_{2}^{b}}{t}$
$\displaystyle=k_{S_{2}}^{+}\Big{\\{}\underbrace{\tfrac{N_{2}}{N_{A}}(P^{b,a}+P^{se,a})}_{\begin{subarray}{c}\text{Total
concentration}\\\ \text{of binding
sites}\end{subarray}}-\underbrace{(S_{2}^{b}+E_{2}^{b}+C_{1}+C_{2})}_{\begin{subarray}{c}\text{Occupied}\\\
\text{binding sites}\end{subarray}}\Big{\\}}S_{2}-k_{S_{2}}^{-}S_{2}^{b}$ (32)
$\displaystyle\;\;\;\;-k_{C_{1}}^{+}S_{2}^{b}E_{1}^{b}+k_{C_{1}}^{-}C_{1},$
$\displaystyle\partialderivative{E_{1}^{b}}{t}$
$\displaystyle=k_{E_{1}}^{+}\Big{\\{}\tfrac{N_{1}}{N_{A}}(P^{b,a}+P^{se,a})-(S_{1}^{b}+E_{1}^{b}+C_{1}+C_{2})\Big{\\}}E_{1}-k_{E_{1}}^{-}E_{1}^{b}$
(33)
$\displaystyle\;\;\;\;-k_{C_{1}}^{+}S_{2}^{b}E_{1}^{b}+(k_{C_{1}}^{-}+k_{C_{1}}^{\text{cat}})C_{1}+\underbrace{k_{C_{2}}^{\text{cat}}C_{2}}_{\begin{subarray}{c}\text{Positive}\\\
\text{feedback}\end{subarray}},$
$\displaystyle\partialderivative{E_{2}^{b}}{t}$
$\displaystyle=k_{E_{2}}^{+}\Big{\\{}\tfrac{N_{2}}{N_{A}}(P^{b,a}+P^{se,a})-(S_{2}^{b}+E_{2}^{b}+C_{1}+C_{2})\Big{\\}}E_{2}-k_{E_{2}}^{-}E_{2}^{b}$
(34)
$\displaystyle\;\;\;\;-k_{C_{2}}^{+}S_{1}^{b}E_{2}^{b}+(k_{C_{2}}^{-}+k_{C_{2}}^{\text{cat}})C_{2}+\underbrace{k_{C_{1}}^{\text{cat}}C_{1}}_{\begin{subarray}{c}\text{Positive}\\\
\text{feedback}\end{subarray}},$ $\displaystyle\partialderivative{C_{1}}{t}$
$\displaystyle=k_{C_{1}}^{+}S_{2}^{b}E_{1}^{b}-(k_{C_{1}}^{-}+k_{C_{1}}^{\text{cat}})C_{1},$
(35) $\displaystyle\partialderivative{C_{2}}{t}$
$\displaystyle=k_{C_{2}}^{+}S_{1}^{b}E_{2}^{b}-(k_{C_{2}}^{-}+k_{C_{2}}^{\text{cat}})C_{2}.$
(36)

### Subendothelium-bound:

$\displaystyle\partialderivative{E_{0}}{t}$
$\displaystyle=-k_{C_{0}}^{+}S_{1}E_{0}+(k_{C_{0}}^{-}+k_{C_{0}}^{\text{cat}})C_{0},\text{
on }\partial\Omega_{\text{inj}},$ (37)
$\displaystyle\partialderivative{C_{0}}{t}$
$\displaystyle=k_{C_{0}}^{+}S_{1}E_{0}-(k_{C_{0}}^{-}+k_{C_{0}}^{\text{cat}})C_{0},\text{
on }\partial\Omega_{\text{inj}}.$ (38)

##  Appendix C Boundary conditions

The boundary conditions for most variables remain consistent across different
geometries. The variables that experience changes in their boundary conditions
when the geometry is altered are those introduced at the domain inlet, such as
the fluid velocity $\overrightarrow{\mathbf{u}}$, the mobile unactivated
platelets $P^{m,u}$, and the kinematic pressure,
$\tilde{\mathrm{p}}=\mathrm{p}/\rho$, at the outlet. The specific conditions
for these variables are detailed in the corresponding subsections specific to
each domain. For all other variables, the following boundary conditions are
applied:

Table 2: Boundary conditions for each variable in the simulation. The table specifies the mathematical boundary condition and the corresponding implementation in OpenFOAM for each boundary. Some variables with the same boundary conditions have been grouped together for brevity.  Variable(s) | inlet | outlet | fixedWalls | injuryWalls
---|---|---|---|---
$\overrightarrow{\mathbf{u}}$ | see below |  | $\partialderivative{\overrightarrow{\mathbf{u}}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\overrightarrow{\mathbf{u}}=0$
---
noSlip
| $\overrightarrow{\mathbf{u}}=0$
---
noSlip
$\tilde{\mathrm{p}}$ |  | $\partialderivative{\tilde{\mathrm{p}}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\tilde{\mathrm{p}}=\mathrm{p}_{\text{out}}$
---
fixedValue
| $\partialderivative{\tilde{\mathrm{p}}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{\tilde{\mathrm{p}}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
$P^{m,u}$ |  | see below
---
codedFixedValue
| $\partialderivative{P^{m,u}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{P^{m,u}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{P^{m,u}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $P^{m,a}$
---
$P^{b,a}$
$P^{b,se}$
| $\partialderivative{P^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{P^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{P^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{P^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
$S_{1}$ |  | $S_{1}=\binom{\text{normal}}{\text{concentration}}$
---
fixedValue
| $\partialderivative{S_{1}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{S_{1}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
|
$-D\partialderivative{S_{1}}{\hat{\mathbf{n}}}=-k_{C_{0}}^{+}S_{1}+k_{C_{0}}^{-}C_{0}$
---
codedMixed
$E_{1}$ |  | $\partialderivative{E_{1}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{E_{1}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{E_{1}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $-D\partialderivative{E_{1}}{\hat{\mathbf{n}}}=k_{C_{0}}^{\text{cat}}C_{0}$
---
codedMixed
$S_{2}$ |  | $S_{2}=\binom{\text{normal}}{\text{concentration}}$
---
fixedValue
| $\partialderivative{S_{2}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{S_{2}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{S_{2}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
$E_{2}$ |  | $\partialderivative{E_{2}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{E_{2}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{E_{2}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{E_{2}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $S_{1}^{b}$
---
$S_{2}^{b}$
$E_{1}^{b}$
$E_{2}^{b}$
$C_{1}$
$C_{2}$
| $\partialderivative{C^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{C^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{C^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{C^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $E_{0}$
---
$C_{0}$
| $\partialderivative{C^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{C^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{C^{k}}{\hat{\mathbf{n}}}=0$
---
zeroGradient
|  See equations (37) and (38)
---
solved in odeSolver.H
ADP |  | $\partialderivative{[\text{ADP}]}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{[\text{ADP}]}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{[\text{ADP}]}{\hat{\mathbf{n}}}=0$
---
zeroGradient
| $\partialderivative{[\text{ADP}]}{\hat{\mathbf{n}}}=0$
---
zeroGradient

###  C.1 Boundary conditions for thrombosis in a rectangular channel

The velocity profile for $\overrightarrow{\mathbf{u}}=[u\;v]^{T}$ is
prescribed using the following function, with a shear rate of
$\dot{\gamma}=1000$ s-1 and radius $r=30\;\mu$m:

$u(y)=-\frac{\dot{\gamma}}{2r}(y-r)^{2}+\frac{1}{2}\dot{\gamma}r.$ (39)

This parabolic profile can be implemented in OpenFOAM using the
codedFixedValue boundary condition

The mobile-unactivated platelets $P^{m,u}$ exhibit margination behavior at the
inlet, where the concentration near the wall is higher compared to the center
of the vessel. To capture this phenomenon, we employ a shape function based on
the inlet profile proposed by Eckstein and Belgacem [24]. The inlet profile
for mobile-unactivated platelets entering a 2D computational domain
$\Omega=[x_{0},\;x_{\text{max}}]\times[y_{0},\;y_{\text{max}}]$ is described
by:

$\displaystyle P^{m,u}(x_{0},y,t)=P_{0}\,c(y),$ (40)

where $P_{0}$ is the normal density of platelets. The shape function $c(y)$ is
defined as:

$\displaystyle c(y)=C_{0}\Big{[}1+KR(y,r)^{m-1}(1-R(y,r))^{n-1}\Big{]}.$ (41)

with $C_{0}$ being a normalizing parameter, $K$ determining the relative
amplitude of the shape, and $m,\;n\in\mathbb{N}$. The function $R(y,r)$ is
given by:

$\displaystyle R(y,r)=\frac{|y-r|}{r},$ (42)

where the vessel radius is $r=(y_{\text{max}}-y_{0})/2$. The exponents are set
as $m=19$ and $n=2$ to produce the effect of a near-wall peak added to a
uniform density $P_{0}$. The normalization parameter is defined as:

$\displaystyle\frac{1}{C_{0}}=\frac{1}{2r}\int_{y_{0}}^{y_{\text{max}}}\Big{[}1+KR(y,r)^{m-1}(1-R(y,r))^{n-1}\Big{]}\,dy$
(43)

Figure 7 shows the specific inlet profile used in the thrombosis example in
Section 4.1, where $r=30\;\mu$m and $K=330$.

Figure 7: Inlet profile for the mobile-unactivated platelets in a vessel with
$r=30\;\mu$m. The peak-to-centre ratios is set to approximately 7.5 using an
amplitude parameter of $K=$ 330\.

###  C.2 Boundary conditions for hemostasis in a microfluidic device

The velocity profile $\overrightarrow{\mathbf{u}}=[u\;v\;w]^{T}$ is prescribed
using the flowRateInletVelocity boundary condition, with volumetric flow rates
$Q_{\text{left}}$ and $Q_{\text{right}}$ at the upper left and upper right
inlets, respectively. The kinematic pressure at the lower outlets is denoted
by $\tilde{\mathrm{p}}_{\text{left}}$ and $\tilde{\mathrm{p}}_{\text{right}}$.
The specific values for the flow rates and kinematic pressure, adapted from
Schoeman et al. [14], are provided in the table below:

Table 3: Fluid velocity and kinematic pressure boundary conditions. Left inlet ($Q_{\text{left}}$) | $1.17\times 10^{-2}\;\text{mm}^{3}/\text{s}$
---|---
Right inlet ($Q_{\text{right}}$) | $9.32\times 10^{-2}\;\text{mm}^{3}/\text{s}$
Left outlet ($\tilde{\mathrm{p}}_{\text{left}}$) | $1.11\;\times 10^{6}\;\text{mm}^{2}/\text{s}^{2}$
Right outlet ($\tilde{\mathrm{p}}_{\text{right}}$) | $1.78\;\times 10^{6}\;\text{mm}^{2}/\text{s}^{2}$

For the mobile-unactivated platelets, the shape function in equation (40) is
adapted for a 3D channel,
$\Omega=[x_{0},\;x_{\text{max}}]\times[y_{0},\;y_{\text{max}}]\times[z_{0},\;z_{\text{max}}]$.
The normalization constant $C_{0}$ is then calculated in two dimensions. The
shape function $c(x,z)$ is defined by:

$\displaystyle
c(x,z)=C_{0}\Big{[}1+K_{x}R(x,r_{x})^{m-1}(1-R(x,r_{x}))^{n-1}+K_{z}R(z,r_{z})^{m-1}(1-R(z,r_{z}))^{n-1}\Big{]}$
(44)

In the 3D hemostasis example, the amplitude parameters are set to
$K_{x}=K_{z}=547$ to create a peak-to-center ratio of approximately 11.9. The
exponents are set to $m=19$ and $n=2$.

##  Appendix D Discretization Schemes

The discretization schemes employed in this work are summarized in Table 4. To
limit spurious oscillations caused by steep gradients at the thrombus edges,
the van Leer flux limiter [25] is applied in the advective terms for platelets
and biochemical species. In Section 3.2, we provided a detailed discussion on
the choice of interpolation scheme for the hindered transport function
$W(\theta^{T})$, which plays a crucial role in the simulation. In the table
below, Theta_Tfa and Theta_Tfd represent the face values of $\theta^{T}$ used
for hindered advection and hindered diffusion, respectively. Both simulations
employ a variable time step $\Delta t$ that is set to ensure the Courant
number, maxCo, does not exceed 0.75. For more practical examples and further
information, please refer to the tutorials directory within the _clotFoam_
repository on GitHub.

Table 4: Discretization schemes for orthogonal meshes in _clotFoam_. ddtSchemes |  |  default |  CrankNicolson 0.9;
---|---
ddt(phi) |  CrankNicolson 0.9;
gradSchemes |  |  default |  Gauss linear;
---|---
grad(p) |  Gauss linear;
divSchemes |  |  default |  none;
---|---
div(phi,U) |  Gauss linear;
div(phiPlt,Plt) |  Gauss limitedVanLeer 0.0 $Pmax;
div(phi,chems) |  Gauss vanLeer;
laplacianSchemes |  |  default |  Gauss linear orthogonal;
---|---
interpolationSchemes |  |  default |  linear
---|---
interpolate(Theta_Tfa) |  downwind phi;
interpolate(Theta_Tfd) |  localMax;
snGradSchemes |  |  default |  orthogonal;
---|---

##  Appendix E Parameter Values

All of the parameters are from various sources as described in Leiderman &
Fogelson [3].

Table 5: Fluid equation parameters. Fluid density ($\rho$) |  1.0 $\times 10^{-3}$ g/mm3
---|---
Dynamic viscosity $(\mu)$ |  2.62507 $\times 10^{-3}$ g/mm/s
Kinematic viscosity ($\nu$) |  2.62507 mm2/s
Carman-Kozeny constant ($C_{CK}$) |  1.0 $\times 10^{6}$ mm-2
Table 6: Diffusion coefficients by species. Platelets ($D_{P}$) | $2.5\times 10^{-5}\;\text{mm}^{2}/\text{s}$
---|---
ADP ($D_{\text{ADP}}$) | $5\times 10^{-4}\;\text{mm}^{2}/\text{s}$
All other chemical species ($D$) | $5\times 10^{-5}\;\text{mm}^{2}/\text{s}$
Table 7: Platelet equation parameters. Platelet diameter $(P_{\text{diam}})$ |  3.0 $\times 10^{-3}\text{ mm}$
---|---
Maximum packing density $(P_{\text{max}})$ |  6.67 $\times 10^{7}$ platelets$/\text{mm}^{3}$
Normal density of platelets $(P_{0})$ |  2.5 $\times 10^{5}$ platelets$/\text{mm}^{3}$
Adhesion rate $(k_{\text{adh}})$ |  3.3212 $\times 10^{-8}\text{ mm}^{3}/\text{s}$
Cohesion rate $(k_{\text{coh}}\times P_{\text{max}})$ |  1.0 $\times 10^{4}\;\text{s}^{-1}$
Threshold for binding affinity $(\eta_{t})$ |  1.0 $\times 10^{-1}$
Rapid change in binding affinity $(\eta^{\ast}+\eta_{t})$ | $0.5-\eta_{t}$
Length of diffusion for the virtual substance $(L_{\eta})$ | $2.0\times P_{\text{diam}}$
Rate of activation by ADP $(k^{\text{pla}}_{\text{adp}})$ |  0.34 $\text{s}^{-1}$
Rate of activation by thrombin $(k^{\text{pla}}_{e_{2}})$ |  0.5 $\text{s}^{-1}$
Critical concentration of ADP $([\text{ADP}]^{\ast})$ |  2.0 $\times 10^{-3}\;\text{nmol of ADP per mm}^{3}$
Critical concentration of thrombin $(e_{2}^{\ast})$ |  1.0 $\times 10^{-6}\;\text{nmol of $e_{2}$ per mm}^{3}$
Total ADP released $(\hat{A})$ |  $2.0\times 10^{-8}$ nmol of ADP per platelet
Total time of ADP secretion ($\tau_{F}$) | 6 s
Step size for calculating ADP secretion ($\Delta\tau$) | 0.25 s
Table 8: Normal concentrations and surface binding site numbers. Chemical | Concentration | Concentration | Binding Sites
---|---|---|---
| (clotFOAM) |  |
$E_{0}$ (TF:VIIa) |  $1.5\times 10^{-7}$ nmol/mm2 |  15 fmol/cm2 |
$E_{1}$ (Xa) |  |  |  $N_{1}$ = 2700
$E_{2}$ (Thrombin) |  |  |  $N_{2}$ = 2000
$S_{1}$ (X) |  $1.7\times 10^{-4}$ nmol/mm3 |  0.17 $\mu$M |  $N_{1}$ = 2700
$S_{2}$ (II) |  $1.4\times 10^{-3}$ nmol/mm3 |  1.4 $\mu$M | $N_{2}=2000$
Table 9: Reactions on subendothelium and platelet surface. Reaction | Complex | Product |  mm3nmol-1s-1 |  M-1 s-1 |  s-1 |  s-1
---|---|---|---|---|---|---
Activation (of-, by-) |  |  | (clotFOAM) |  |  |
($S_{1}$, $E_{0}$) | $C_{0}$ | $E_{1}$ | $k_{C_{0}}^{+}=8.95\times 10^{3}$ | $k_{C_{0}}^{+}=8.95\times 10^{6}$ | $k_{C_{0}}^{-}=1.0$ | $k_{C_{0}}^{\text{cat}}=1.15$
($S_{2}^{b}$, $E_{1}^{b}$) | $C_{1}$ | $E_{2}^{b}$ | $k_{C_{1}}^{+}=1.03\times 10^{5}$ | $k_{C_{1}}^{+}=1.03\times 10^{8}$ | $k_{C_{1}}^{-}=1.0$ | $k_{C_{1}}^{\text{cat}}=30.0$
($S_{1}^{b}$, $E_{2}^{b}$) | $C_{2}$ | $E_{1}^{b}$ | $k_{C_{2}}^{+}=1.73\times 10^{4}$ | $k_{C_{2}}^{+}=1.73\times 10^{7}$ | $k_{C_{2}}^{-}=1.0$ | $k_{C_{2}}^{\text{cat}}=0.23$
Table 10: Binding on platelet surface Reaction | Reactants | Product |  mm3nmol-1s-1 |  M-1 s-1 |  s-1
---|---|---|---|---|---
|  |  | (clotFOAM) |  |
$S_{1}$ (V) |  $S_{1}$, $P_{1}$ | $S_{1}^{b}$ | $k_{S_{1}}^{+}=5.7\times 10^{4}$ | $k_{S_{1}}^{+}=5.7\times 10^{7}$ | $k_{S_{1}}^{-}=0.17$
$E_{1}$ (Xa) |  $E_{1}$, $P_{1}$ | $E_{1}^{b}$ | $k_{E_{1}}^{+}=1.0\times 10^{4}$ | $k_{E_{1}}^{+}=1.0\times 10^{7}$ | $k_{E_{1}}^{-}=2.5\times 10^{-2}$
$S_{2}$ (II) |  $S_{2}$, $P_{2}$ | $S_{2}^{b}$ | $k_{S_{2}}^{+}=1.0\times 10^{4}$ | $k_{S_{2}}^{+}=1.0\times 10^{7}$ | $k_{S_{2}}^{-}=5.9$
$E_{2}$ (IIa) |  $E_{2}$, $P_{2}$ | $E_{2}^{b}$ | $k_{E_{2}}^{+}=1.0\times 10^{4}$ | $k_{E_{2}}^{+}=1.0\times 10^{7}$ | $k_{E_{2}}^{-}=5.9$

Supporting Information

## Mesh generation for 3D hemostasis example

To define the mesh in OpenFOAM, we utilize 3 blocks as seen in Figure 8. The
creation of holes in the vertical channels at the interface with the injury
channel is accomplished using the mergePatchPairs function. The dimensions of
the H-domain, which serve as the computational domain, are provided in Table
11, following the description given in Schoeman et al. [14]. The mesh is
defined so that the cells in the injury channel have a uniform spacing $\Delta
x=\Delta y=\Delta z=$2.5 $\mu$m. The mesh outside of the injury channel is
graded for computational efficiency. The resulting mesh, as shown in Figure 8,
consists of 37,248 finite volume cells. The specific grading parameters used
to generate the mesh can be found in the file
tutorials/Hjunction3D/system/blockMeshDict in the _clotFoam_ repository on
GitHub. These parameters were calculated using the free blockMesh grading tool
[26].

Table 11: Dimensions of the computational domain of the H-shaped microfluidic device. Width of vertical channel ($w_{\text{chan}}$) |  100 $\mu$m
---|---
Height of vertical channel ($h_{\text{chan}}$) |  50 $\mu$m
Length of vertical channel ($\ell_{\text{chan}}$) |  150 $\mu$m
Width of injury channel ($w_{\text{inj}}$) |  50 $\mu$m
Height of injury channel ($h_{\text{inj}}$) |  20 $\mu$m
Length of injury channel ($\ell_{\text{inj}}$) |  150 $\mu$m

Figure 8: Meshing the H-shaped microfluid device.

## Mesh convergence study

To evaluate the accuracy and convergence of the _clotFoam_ solver, we
performed a mesh convergence study in a rectangular domain measuring 120
$\mu$m in length, 30 $\mu$m in height, and with an injury length of 45 $\mu$m.
The objective of this study was to assess the error by comparing the solution
between a fine mesh and a coarse mesh. We utilized five different mesh
discretizations, with uniform cell heights and widths ($h$) of 3, 1.5, 0.75,
0.375, and 0.1875 $\mu$m. In this study, the finest mesh was denoted as
$h_{\text{fine}}=0.1875\;\mu$m, and its corresponding solution was represented
by $\psi_{\text{fine}}$. To facilitate the comparison of solutions across
different mesh resolutions, we employed the Resample with Dataset filter in
ParaView to interpolate the coarse meshes onto the finest mesh. To ensure
numerical stability, the time step was selected as a function of $h$ in a way
that the maximum Courant number ($C_{0}$) did not exceed 0.75 for each mesh.
The Courant number is defined as $C_{0}=\frac{\Delta
t\|\overrightarrow{\mathbf{u}}\|}{h}$, where $\Delta t$ represents the time
step and $\|\overrightarrow{\mathbf{u}}\|$ is the magnitude of the velocity
field.

The relative error for each solution $\psi_{h_{k}}$ with mesh $h_{k}$ is
calculated using the following formula:

$E(\psi_{h_{k}})=\frac{\|\psi_{\text{fine}}-\psi_{h_{k}}\|_{\ell_{2}}}{\|\psi_{\text{fine}}\|_{\ell_{2}}},$
(45)

where the $\ell_{2}$ norm is defined for a solution vector
$\overrightarrow{\mathbf{x}}\in\mathbb{R}^{n}$ as:

$\displaystyle\|\overrightarrow{\mathbf{x}}\|_{\ell_{2}}$
$\displaystyle=\sqrt{\sum_{i=1}^{n}x_{i}^{2}},$ (46)

The convergence rate, denoted by $r$, was approximated by fitting a line to
the error vectors on a log-log scale using least squares.

In Figure 9, we present the errors at time $t=5$ seconds for four different
species that represent various components of the clotting model:

  1. 1.

$S_{1}$ is a fluid-phase biochemical species that enters the domain with a
fixed concentration and reacts on subendothelium with $E_{0}$ to produce the
enzyme $E_{1}$.

  2. 2.

$E_{1}$ is a fluid-phase species that enters the domain through a reactive
boundary condition on the subendothelium.

  3. 3.

$P^{m,u}$ is fluid-phase platelet species that enters the domain with a
marginated profile and experiences hindered transport in the presence of a
growing thrombus. This species has reaction terms that depend on the adhesion
region and biochemical agonists such as ADP and thrombin.

  4. 4.

$E_{2}+E_{2}^{b}$ represents the total thrombin concentration and serves as a
surrogate measure of overall clot formation. $E_{2}$ is a fluid-phase species,
while $E_{2}^{b}$ is platelet-bound.

The error plot in Figure 9 demonstrates that the fluid-phase biochemical
species converge with second-order accuracy, while the platelet-related
species exhibit a super-linear convergence rate.

Figure 9: Mesh convergence study in a rectangular domain. Errors are
calculated at $t=5$ seconds with respect to a fine mesh solution
$\psi_{\text{fine}}$ ($h=0.1875\;\mu$m), which serves as a reference solution.
The approximate convergence rate $r$ is displayed in the legend for each
species.

## Effect of mesh topology

To investigate the impact of mesh topology on thrombus formation, we generated
an unstructured mesh using the open-source software Gmsh, as shown in Figure
10. Our preliminary findings suggest that the specific mesh topology has a
limited effect on the overall growth of the thrombus. However, we observed
variations in clot density, with less dense regions occurring in areas where
the mesh topology is more skewed. It is important to note that these
observations are based on a single test with an unconventional mesh, and
further investigations incorporating carefully selected discretization schemes
may help mitigate this issue.

To assess the influence of mesh topology on thrombus growth, we conducted a
comparative study using two different mesh configurations within a rectangular
domain, as depicted in Figure 10. The left column represents simulations
conducted using a structured mesh, while the right column corresponds to
simulations conducted using an unstructured mesh. Our comparison revealed
similar patterns of thrombus growth in both simulations. We quantitatively
evaluated the differences in thrombi by calculating the integral of the bound
platelet fraction over the domain and normalizing it by its maximum value,
resulting in the “normalized amount of bound platelets.” The bottom plots in
Figure 10 demonstrate that the error in the normalized amount of bound
platelets is less than 7% between the structured and unstructured mesh
configurations after 600 seconds of simulation.

Figure 10: Comparison of the effect of mesh topology on thrombus growth for a
rectangular domain measuring 240 $\mu$m by 60 $\mu$m. The top row depicts a
structured (uniform) mesh on the left and an unstructured mesh on the right.
The second row provides a close-up view of the bound platelet fraction in a
120 $\mu$m long by 30 $\mu$m high region surrounding the thrombus after 600 s
of simulation. The bottom row quantifies the discrepancy between the
simulations by integrating the bound platelet fraction over the entire domain
and normalizing it by the maximum value of that quantity.

## Quantitative validation

To quantitatively validate the _clotFoam_ solver, we conducted comparative
studies with the works of Leiderman and Fogelson [3] and Schoeman et al. [14].
It is important to acknowledge that the clotting model implemented in
_clotFoam_ is a significantly reduced version of the Leiderman and Fogelson
model. This reduction was intended to provide researchers with a flexible
framework that can be adapted to suit their specific research needs and
objectives. As a result, direct comparisons between these models can be
challenging due to inherent differences in their formulations. However,
considering the results obtained from these comparative studies, along with
the mesh convergence analysis depicted in Figure 9, it becomes evident that
_clotFoam_ demonstrates reliability in simulating thrombus growth under
arterial flow conditions.

#### Comparison with Leiderman and Fogelson 2011

In the first comparative study, we compared thrombus growth from the
_clotFoam_ solver to Figure 8 in Leiderman and Fogelson [3]. We simulated
thrombus growth under three shear rates (500, 1000, and 1500 s-1) with four
different inlet platelet profiles ($P_{0}$, $P_{1}$, $P_{2}$, $P_{3}$ as
defined in Figure 3 of Leiderman and Fogelson) that have increasing
margination ratios. To quantify the differences, we computed the ’area’ in
which the bound platelet concentration exceeded 10%, 50%, and 90% of
$P_{\text{max}}$, respectively. The first percentage provides a measure of the
overall thrombus size, while the latter two percentages provide information
about the bound platelet density distribution within the thrombus. These
’areas’ were computed in ParaView by counting the grid cells in which the
various platelet concentration levels were exceeded, and are presented in
Figure 11.

Figure 11: Quantitative analysis of clot ’area’ for comparison with Fig. 8 in
Leiderman & Fogelson [3]. The top row (a-c) presents results obtained using a
uniform platelet profile. The middle row (d-f) shows results obtained using a
fixed nonuniform platelet profile $P_{1}$, as defined in Fig. 3 of Leiderman &
Fogelson [3]. The bottom row (g–i) displays results obtained using platelet
profiles $P_{1}$, $P_{2}$, and $P_{3}$ for wall shear rates of 500, 1000, and
1500 s-1, respectively. In all cases, the dotted line, dashed line, and solid
line represent wall shear rates of 500, 1000, and 1500 s-1, respectively.

Similar to the plots in Leiderman and Fogelson [3], we observed the appearance
of a temporary plateau at an area of 100 grid cells. This corresponds to the
region where the adhesion rate function $k_{\text{adh}}H_{\text{adh}}(x)$ is
non-zero. However, for the highest concentration of bound platelets, when the
bound platelets concentration exceeded 90% of $P_{\text{max}}$, we observed a
plateau at an area of 50 grid cells. This indicates that only half of the
adhesion region exceeds a density of 90% of $P_{\text{max}}$. This behavior
can be seen in Figure 5 of our manuscript, where the adhesion region at the
bottom of the thrombus has a less dense concentration compared to the rest of
the clot. Given the reduced nature of the _clotFoam_ model, it is not
surprising that our results for the timing of plateaus, overall clot growth,
and shear rate dynamics differ from those published by Leiderman and Fogelson.
However, we do observe similar characteristics of thrombus growth, with the
curves in Figure 11 following the same general trajectories as those reported
by Leiderman and Fogelson.

#### Comparsion with Schoeman et al. 2016

In the second comparative study, we conducted an assessment of thrombus
formation over time by comparing our results with Figure 5 of Schoeman et al.
[14]. Our primary objective with the illustrative example of hemostasis in a
microfluidic device using the _clotFoam_ software was to demonstrate its
capability in simulating clot formation within such a device. However, it is
important to acknowledge that the platelet model employed by _clotFoam_ does
not incorporate crucial factors like von Willebrand factor (vWF) or shear
dependence in the platelet aggregation process, which are known to play
significant roles in high shear settings, such as those encountered in the
H-shaped microfluidic device described by Schoeman et al.

Despite the limitations of the platelet aggregation model in _clotFoam_ ,
Figure 12 illustrates characteristic platelet accumulation in the injury
channel for whole blood on collagen and tissue factor (TF). The plot shows
that as time progresses, the platelet accumulation increases until the injury
channel reaches full occlusion at approximately 30 minutes, resulting in a
plateau in the amount of bound platelets. This behavior is similar to what is
observed in the Schoeman et al. paper, although the time scales differ. These
findings demonstrate that the base clotting model in _clotFoam_ is capable of
simulating hemostasis in a microfluidic device. Researchers can easily modify
the model to incorporate the effects of shear, von Willebrand factor (vWF),
thromboxane A2, fibrin, and other relevant biochemical species involved in the
coagulation process. The primary objective of developing _clotFoam_ was to
provide researchers with a versatile framework that can be easily customized
and adapted to different clotting models. This flexibility allows for further
exploration and customization to meet the specific research requirements of
various studies in the field of hemostasis and thrombosis.

Figure 12: Characteristic platelet accumulation in the injury channel for
whole blood on collagen-TF. The plot shows the progressive accumulation of
platelets over time until the injury channel reaches occlusion at
approximately 30 minutes.

## References

  * [1] A. L. Fogelson, A. L. Kuharsky, Membrane binding-site density can modulate activation thresholds in enzyme systems, Journal of Theoretical Biology 193 (1) (1998) 1–18.  doi:10.1006/jtbi.1998.0670.
  * [2] A. L. Kuharsky, A. L. Fogelson, Surface-Mediated Control of Blood Coagulation: The Role of Binding Site Densities and Platelet Deposition, Biophysical Journal 80 (3) (2001) 1050–1074.  doi:10.1016/S0006-3495(01)76085-7.
  * [3] K. Leiderman, A. L. Fogelson, Grow with the flow: a spatial–temporal model of platelet deposition and blood coagulation under flow, Mathematical medicine and biology: a journal of the IMA 28 (1) (2011) 47–84.  doi:10.1093/imammb/dqq005.
  * [4] K. Miyazawa, A. L. Fogelson, K. Leiderman, Inhibition of platelet-surface-bound proteins during coagulation under flow i: Antithrombin and heparin, Biophysical journal 122 (1) (2023) 230–240.  doi:10.1016/j.bpj.2022.11.023.
  * [5] K. B. Neeves, K. Leiderman, Mathematical models of hemostasis, Trauma induced coagulopathy (2016) 567–584doi:10.1007/978-3-319-28308-1\\_35.
  * [6] K. Leiderman, B. Bannish, M. Kelley, A. Palmisano, Mathematical models of thrombus formation and fibrinolysis, in: Cardiovascular thrombus: from pathology and clinical presentation to imaging, pharmacotherapy and interventions, Academic Press San Diego, 2018, pp. 67–82.  doi:10.1016/B978-0-12-812615-8.00005-3.
  * [7] S. L. Diamond, Systems biology of coagulation, Journal of Thrombosis and Haemostasis 11 (2013) 224–232.  doi:10.1111/jth.12220.
  * [8] S. Yesudasan, R. D. Averett, Recent advances in computational modeling of fibrin clot formation: A review, Computational biology and chemistry 83 (2019) 107148.  doi:10.1016/j.compbiolchem.2019.107148.
  * [9] M. Anand, M. Panteleev, F. Ataullakhanov, Computational models of hemostasis: Degrees of complexity, Applications in Engineering Science 10 (2022) 100103.  doi:10.1016/j.apples.2022.100103.
  * [10] J. O. Taylor, R. S. Meyer, S. Deutsch, K. B. Manning, Development of a computational model for macroscopic predictions of device-induced thrombosis, Biomechanics and modeling in mechanobiology 15 (2016) 1713–1731.  doi:10.1007/s10237-016-0793-2.
  * [11] V. Govindarajan, S. Zhu, R. Li, Y. Lu, S. L. Diamond, J. Reifman, A. Y. Mitrophanov, Impact of tissue factor localization on blood clot structure and resistance under venous shear, Biophysical journal 114 (4) (2018) 978–991.  doi:10.1016/j.bpj.2017.12.034.
  * [12] R. Méndez Rojano, M. Zhussupbekov, J. F. Antaki, D. Lucor, Uncertainty quantification of a thrombosis model considering the clotting assay pfa-100®, International Journal for Numerical Methods in Biomedical Engineering 38 (5) (2022) e3595.  doi:doi.org/10.1002/cnm.3595.
  * [13] A. Bouchnita, A. V. Belyaev, V. Volpert, Multiphase continuum modeling of thrombosis in aneurysms and recirculation zones, Physics of Fluids 33 (9) (2021) 093314.  doi:10.1063/5.0057393.
  * [14] R. M. Schoeman, K. Rana, N. Danes, M. Lehmann, J. A. Di Paola, A. L. Fogelson, K. Leiderman, K. B. Neeves, A microfluidic model of hemostasis sensitive to platelet function and coagulation, Cellular and molecular bioengineering 10 (2017) 3–15.  doi:10.1007/s12195-016-0469-0.
  * [15] N. A. Danes, K. Leiderman, A density-dependent fem-fct algorithm with application to modeling platelet aggregation, International journal for numerical methods in biomedical engineering 35 (9) (2019) e3212.  doi:10.1002/cnm.3212.
  * [16] K. Leiderman, A. L. Fogelson, The influence of hindered transport on the development of platelet thrombi under flow, Bulletin of mathematical biology 75 (2013) 1255–1283.  doi:10.1007/s11538-012-9784-3.
  * [17] M. Rezaeimoghaddam, F. N. van de Vosse, Continuum modeling of thrombus formation and growth under different shear rates, Journal of Biomechanics 132 (2022) 110915.  doi:10.1016/j.jbiomech.2021.110915.
  * [18] W.-T. Wu, M. A. Jamiolkowski, W. R. Wagner, N. Aubry, M. Massoudi, J. F. Antaki, Multi-constituent simulation of thrombus deposition, Scientific reports 7 (1) (2017) 1–16.  doi:10.1038/srep42720.
  * [19] C. Greenshields, OpenFOAM v9 User Guide, The OpenFOAM Foundation, London, UK, 2021.
  * [20] K. N. Shankar, Y. Zhang, T. Sinno, S. L. Diamond, A three-dimensional multiscale model for the prediction of thrombus growth under flow with single-platelet resolution, PLOS Computational Biology 18 (1) (2022) e1009850.  doi:10.1371/journal.pcbi.1009850.
  * [21] R. Méndez Rojano, A. Lai, M. Zhussupbekov, G. W. Burgreen, K. Cook, J. F. Antaki, A fibrin enhanced thrombosis model for medical devices operating at low shear regimes or large surface areas, PLOS Computational Biology 18 (10) (2022) e1010277.  doi:10.1371/journal.pcbi.1010277.
  * [22] R. I. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, Journal of computational physics 62 (1) (1986) 40–65.  doi:10.1016/0021-9991(86)90099-9.
  * [23] M. Félétou, The endothelium, part i: Multiple functions of the endothelial cells – focus on endothelium-derived vasoactive mediators, Colloquium Series on Integrated Systems Physiology: From Molecule to Function 3 (4) (2011) 1–306.  doi:10.4199/c00031ed1v01y201105isp019.
  * [24] E. C. Eckstein, F. Belgacem, Model of platelet transport in flowing blood with drift and diffusion terms, Biophysical journal 60 (1) (1991) 53–69.  doi:10.1016/S0006-3495(91)82030-6.
  * [25] B. van Leer, Towards the ultimate conservative difference scheme. ii. monotonicity and conservation combined in a second-order scheme, Journal of Computational Physics 14 (4) (1974) 361–370.  doi:10.1016/0021-9991(74)90019-9.
  * [26] Scripts/blockmesh grading calculation, https://openfoamwiki.net/index.php/Scripts/blockMesh_grading_calculation (2020).