A simplified mesoscale 3D model for characterizing fibrinolysis under flow conditions

One of the routine clinical treatments to eliminate ischemic stroke thrombi is injecting a biochemical product into the patient’s bloodstream, which breaks down the thrombi’s fibrin fibers: intravenous or intravascular thrombolysis. However, this procedure is not without risk for the patient; the worst circumstances can cause a brain hemorrhage or embolism that can be fatal. Improvement in patient management drastically reduced these risks, and patients who benefited from thrombolysis soon after the onset of the stroke have a significantly better 3-month prognosis, but treatment success is highly variable. The causes of this variability remain unclear, and it is likely that some fundamental aspects still require thorough investigations. For that reason, we conducted in vitro flow-driven fibrinolysis experiments to study pure fibrin thrombi breakdown in controlled conditions and observed that the lysis front evolved non-linearly in time. To understand these results, we developed an analytical 1D lysis model in which the thrombus is considered a porous medium. The lytic cascade is reduced to a second-order reaction involving fibrin and a surrogate pro-fibrinolytic agent. The model was able to reproduce the observed lysis evolution under the assumptions of constant fluid velocity and lysis occurring only at the front. For adding complexity, such as clot heterogeneity or complex flow conditions, we propose a 3-dimensional mesoscopic numerical model of blood flow and fibrinolysis, which validates the analytical model’s results. Such a numerical model could help us better understand the spatial evolution of the thrombi breakdown, extract the most relevant physiological parameters to lysis efficiency, and possibly explain the failure of the clinical treatment. These findings suggest that even though real-world fibrinolysis is a complex biological process, a simplified model can recover the main features of lysis evolution.


B Lysis of the n-th slice with the analytical model
Before initiating the lysis of the second slice x ∆ 1 , the quantity F 0 u f ∆ t has to be advected over a distance ∆. Hence, the anti-FA that accumulates in x ∆ 1 is given by: The Heaviside step function Θ reflects that anti-FA does not reach x ∆ 1 until x ∆ 0 has been lysed, and the anti-FA has been transported over the distance ∆. Thus, at the time t = t 0 + ∆ u f , the concentration of anti-FA at x ∆ 1 is F 0

1/10
After injecting (S2) into (1) and integrating, we have : which after rearrangement gives: Keeping the positive solution of the resulting quadratic equation, we get the time t 1 at which the second slice x ∆ 1 is lysed: By applying the same reasoning recursively, we get the formula for the time to lyse n slices :

C Analytical model independent of ∆
We show in Figure S2 that the analytical formula for the lysis time of the thrombus does not depend on the choice of the spatial slicing ∆, provided that ∆ ≪ L. We demonstrate it here analytically, starting with equation (5) and rearranging terms: By taking the first order approximation of ∆, we neglect the second order term ∆ 2 u 2 f , which yields: On the other hand, we can write : 10 -5 10 -4 10 -3 Figure S2. Lysis profiles of a 1 cm thrombus predicted with the analytical formula. F 0 = 3.2 mg·ml −1 , F 0 = 10 ng·ml −1 , ∇P = 10 6 Pa·m −1 , k 1 = 12.5(mg·ml −1 ) −1 . The lysis time predicted with the analytical formula is independent of the slicing of the thrombus ∆, is ∆ is small compared to L.
where dt is the variation of lysis time at time t for the slice ∆. At first order in ∆, t n + t n−1 = 2t is a sufficient approximation. Combining (S5) and (S6), and writing ∆ = dh we get: By integrating, and setting h = 0 at t = 0 we finally get: which, using the properties of the log, can be rewritten as: or: and which is independent of the slicing.

D Experimental lysis fronts with various concentrations.
We put at disposition here lysis front measurements for fibrinolysis experiments that we conducted, while varying initial fibrin concentration ( Figure S3A), and the tPA solution concentration ( Figure S3B). The dots represent the experimental

E Pressure boundary conditions
Pressure boundary conditions (BCs) were imposed on both ends of the tube by imposing densities, as is usually done in LB. The 3D Zou-He pressure BC 1 and a novel pressure BC were tested. The novel BC showed better simulation stability when increasing the pressure difference; that is why it was chosen for the numerical model. This novel BC can also simulate a pulsatile pressure profile by applying a time modulation to the inlet density. Figure S4 presents an explanation of how the populations at the boundaries are computed. Let us show how the density is imposed with this pressure BC. We consider the situation displayed in Figure S4, where we need to impose the density on the outlet boundary node N. After having copied the unknown populations f i from node N − 1, the density at node N can be computed as usually in LB: where the f j designates the known populations from the streaming step, at node N. The equilibrium function f eq k at node N, used in the standard BGK collision, reads 2 : where w k are the weighting factors, c k the lattice discrete velocities, and c s the speed of sound magnitude in the lattice. By using the property that ∑ k f eq k (x N ,t) = ∑ k f k (x N ,t) = ρ tmp at node N, we subtract and add the following equilibrium functions

4/10
to the populations f k at node N: By juggling with these equilibrium functions, we can thus impose the density ρ target at the boundary. It is worth noting that the velocities u N remain unchanged.

F Particles transport
The particles in the numerical model embody the anti-FA entity, which accounts for the lysis of fibrin. They are subject to advection, the transport driven by the fluid (lytic solution), and diffusion, which is the movement due to thermal agitation and molecular collisions.

F.1 Advection
The Palabos library offers a built-in algorithm that automatically transports particles advected by the fluid. However, it is a one-size-fits-all solution and is not necessarily the most efficient way to advect the particle for a specific model. Looping over all particles at each iteration step is costly. In our thrombolysis model, the flow varies very slowly during most of the lysis. It is thus interesting in our case, to increase the advection time scale according to the flow velocity. Practically, instead of moving on a distance dx the particles at each lattice-Boltzmann iteration, we move them of T dx after T iterations, where T is computed such that the particles move at most of half a lattice site: . (S15) This computational trick allows a speedup of around a factor 10 in a typical lysis simulation, without any visible difference on the lysis curves. The advection velocity of each particle is computed as with the Palabos solver, using weighted averaging of the velocity at the 4 closest lattice nodes.

F.2 Diffusion
The Péclet number: where u is the mean velocity of the fluid, L is a characteristic dimension, and D the diffusion coefficient of the diffusing mass in the medium, gives an approximation of the ratio of advection against diffusion. For high Péclet (Pe ≫ 1), diffusion can be neglected but should otherwise be taken into account. In the clotted tubes of the experiments presented in the manuscript, and for tPA molecules, the initial Péclet is: which shows that the diffusion of tPA is non-negligible to initiate the lysis of the clot.
Since the anti-FA is modeled by particles, the most straightforward way to implement their diffusion is by individually giving a random velocity at each simulation time step dt, of norm : where n = 3 is the spatial dimension of the model. It can be verified that this method yields the expected collective behavior, by computing the mean squared displacement (MSD) of the N particles: is the 3D position of particle i at time t.
It can be shown that the MSD approximates the diffusion coefficient: By simulating N = 100 particles in a system with no advection, already good agreement can be observed with eq. (S20), as shown in Figure S5.

G Porosity and permeability validation
Let us here validate our implementation of the PBB model for computing flow in simulated fibrin clots. Let us recall Darcy's law for laminar flow 3 : with ν the kinematic viscosity, ρ the fluid density, ∆P the pressure drop and k the permeability. Walsh et al. show that the permeability k PBB in m 2 of PBB voxels is related to γ and the simulation time step ∆t as follows: Combining (14), (S21) and (S22), we can express Darcy's law in terms of γ: We measured the values of the fluid velocity u f in numerical simulations, and compared with the analytical solution of Darcy and Walsh equation (S23) in Figure S6A. Flow was simulated in a 3D rectangular duct filled with a porous medium with bounce-back fractions of γ = 0.2, 0.4, 0.8 respectively. The velocity profile was measured along the y-axis (transversal) in the center of the duct. In the analytical solution, we imposed no-slip conditions, i.e. u f = 0 at the walls. Next, we investigate the validity of the PBB method for a typical range of fibrin clot permeabilities. Other commonly used laws to describe fibrous porous media permeabilities are Clague's equation 4 , a numerical solution of permeability of randomly organized fibrous media: and Jackson-James' equation 5 , the weighted average of the solution to the Stokes equation for flow parallel to or normal to 2D periodic square arrays of cylinders: Permeabilities measured with the model as a function of input solid fractions are shown in Figure S6B-D, and compared with their empirical counterparts equations (16),(S24) and (S25). The permeability values for the physiological range of low solid fractions (0.003 − 0.04, i.e. 1 − 4 mg·ml −1 fibrinogen concentration) do not vary substantially from one model to another; thus, if the clots studied are highly porous, the choice of model has little impact. Clague's law predicts an increase of permeability with the increase of solid fraction in the 0.8 − 1 range, and Jackson-James' law asymptotically approaches zero as n * s approaches 0.4. Thus, in order to be able to simulate the whole range of solid fractions, Davies' permeability law was chosen for the model.

H Deviation of numerical from analytical model
We show here the deviation of the mesoscopic numerical lysis front, from the front given by the analytical formula, in the conditions of Figure 4.The relative absolute difference ε on the front position in time is shown on Figure S7. We see that the deviation stays below 5% during most of the lysis, and increases up to 20% in the last third of the total lysis time. This increase in the difference is due to the assumption of constant flow made in the analytical model.

I Lysis of realistic heterogeneous clots
We generated numerical thrombi directly from ischemic stroke patient thrombi retrieved by thrombectomy and sliced by Staessens et al. 8 . Platelet stained slices are converted to grayscale using the 'L' algorithm of the Python Imaging Library.  5 , according to equations (16),(S24) and (S25) respectively. Simulations were made in a tube of radius 4.5 mm (23 nodes), length 1 cm (46 nodes), with a fully occluding homogeneous thrombus (i.e., same dimensions as the tube). The space-time discretization was ∆x = 2.17 · 10 −4 m and ∆t = 1.58 · 10 −3 s, and the pressure gradient ∇P = 100 Pa·m −1 The hydrated fiber radius ranges from R f = 51 to 60 nm for formulae and simulations, as in 6 .
We then assign a fibrin concentration in a linear scale ranging from 1 mg·ml −1 (white, no platelets) to 8 mg·ml −1 (black, platelet-rich), Figure S8A. Platelet-rich zones are supposed to be more tPA-resistant and to have lower permeability 9 , 10 , 11 , which we translate in the present model with higher fibrin concentration regions. We lysed a RBC-rich clot and a PLT-rich clot with the numerical model, as well as homogeneous clots with the same average concentration of 2.45 and 3.92 mg·ml −1 respectively, Figure S8B. The lysis of the heterogeneous clots was slightly faster, and the recanalization occurred slightly before, Figure S8C  The mean throughput is monitored during the lysis. RBC-rich clots are more permeable since before the start of the lysis, the average fibrin concentration being lower than the PLT-rich's one.