A Mathematical Model of Platelet Aggregation in an Extravascular Injury Under Flow
Kathryn G. Link, Matthew G. Sorrells, Nicholas A. Danes, Keith B. Neeves, Karin Leiderman, Aaron L. Fogelson
AA MATHEMATICAL MODEL OF PLATELET AGGREGATION INAN EXTRAVASCULAR INJURY UNDER FLOW. ∗ KATHRYN G. LINK † , MATTHEW G. SORRELLS ‡ , NICHOLAS A. DANES § , KEITH B.NEEVES ¶ , KARIN LEIDERMAN § , AND
AARON L. FOGELSON †(cid:107)
Abstract.
We present the first mathematical model of flow-mediated primary hemostasis inan extravascular injury which can track the process from initial deposition to occlusion. The modelconsists of a system of ordinary differential equations (ODE) that describe platelet aggregation(adhesion and cohesion), soluble-agonist-dependent platelet activation, and the flow of blood throughthe injury. The formation of platelet aggregates increases resistance to flow through the injury, whichis modeled using the Stokes-Brinkman equations. Data from analogous experimental (microfluidicflow) and partial differential equation models informed parameter values used in the ODE modeldescription of platelet adhesion, cohesion, and activation. This model predicts injury occlusionunder a range of flow and platelet activation conditions. Simulations testing the effects of shear andactivation rates resulted in delayed occlusion and aggregate heterogeneity. These results validate ourhypothesis that flow-mediated dilution of activating chemical ADP hinders aggregate development.This novel modeling framework can be extended to include more mechanisms of platelet activationas well as the addition of the biochemical reactions of coagulation, resulting in a computationallyefficient high throughput screening tool.
Key words.
Mathematical modeling, flows in porous media, mathematical biology, bloodclotting, hemostasis
AMS subject classifications.
1. Introduction.
Hemostasis is the first line of defense upon vascular injury (ei-ther a disruption in the subendothelial lining or a breach in the vessel wall) wherebya blood clot (thrombus) forms to prevent the loss of blood [47, 48]. The hemostaticresponse consists of two intertwined processes: platelet aggregate formation and co-agulation [13, 34]. These processes are triggered when reactive proteins either inthe extravascular matrix or on the subendothelium are exposed to the blood plasma.Platelets are blood cells that circulate in the human vasculature in their unactivatedstate and mediate the biophysical and biochemical aspects of thrombus formation.Platelets become activated when in contact with collagen in the subendothelium orwith chemical agonists in the plasma. Activated platelet surfaces support coagulationreactions that produce the enzyme thrombin. Thrombin in turn cleaves the solubleblood protein fibrinogen into fibrin monomers, which polymerize to form a gel. Thisgel provides structure and stability to the aggregate. The size and structure of an ag-gregate as well as the time it takes the aggregate to form, depend not only on plateletfunction and coagulation, but on the local hydrodynamic environment. Blood flowingin a vessel is subject to a pressure difference across the vascular wall and any dis-ruption of the wall can lead to loss of blood, which greatly affects the delivery andremoval of platelets, chemical agonists, and coagulation proteins. ∗ Submitted on March 4, 2020.
Funding:
This work was funded by the National Institutes of Health under grants R01HL120728and R61HL141794 and National Science Foundation grant CBET-1351672. † ∼ link). ‡ Department of Chemical and Biological Engineering, Colorado School of Mines, Golden, CO. § Department of Applied Mathematics and Statistics, Colorado School of Mines, Golden, CO. ¶ Departments of Bioengineering and Pediatrics, Hemophilia and Thrombosis Center, Universityof Colorado Denver, Anschutz Medical Campus, Aurora, CO (cid:107)
Department of Biomedical Engineering University of Utah, Salt Lake City, UT.1
This manuscript is for review purposes only. a r X i v : . [ q - b i o . CB ] M a r K.G. LINK, ET AL.
Upon extravascular injury, two important reactive proteins are exposed: collagenand tissue factor (TF). Blood cells and plasma containing coagulation proteins leakfrom the vessel into the extravascular space. Platelets flow into the injury and beginto deposit; they adhere to the collagen, and bind other insoluble agonists, which inturn triggers the activation of key integrins on platelet surfaces [2, 4, 14, 21, 48]. Thesoluble agonist ADP, which is released from platelet granules into the fluid by activatedplatelets themselves, interacts with platelet receptors P2Y and P2Y to additionallytrigger exposure of surface integrins [3, 18, 33]. This effect of soluble agonist ADPis an example of activation without contact with the subendothelial matrix [1, 27].Activated platelets provide surfaces to which more platelets can cohere, enhancingto the growth of the platelet aggregate to prevent further blood leakage. Plateletaggregate formation depends strongly on both the hydrodynamic environment andthe formation and dissociation of bonds that enable platelet adhesion and cohesion.Failure of any of the processes that mediate platelet aggregation, under a variety offlow conditions, can result in significant blood loss.In recent years there has been an increased focus on understanding the molecu-lar basis of bleeding disorders and associated variability in levels of key coagulationfactors, platelet count, and insoluble proteins. Most clinical assays used to assessbleeding risk are performed under static conditions. However, data has shown thathydrodynamic forces affect platelet aggregation and fibrin formation, both key com-ponents of a thrombus [31]. This has motivated the development of mathematicalmodels of hemostasis that integrate platelet function, coagulation, and hemodynam-ics to predict bleeding risk based on measurable biochemical and biophysical factors.Our group has developed mathematical models of flow-mediated coagulation inan intravascular injury setting to predict regulatory mechanisms of hemostasis andthrombosis [10, 12, 15, 22, 23, 24, 25, 28, 29]. These models use a continuum approachand express dynamics using differential equations. Our ordinary differential equation(ODE) models of thrombosis employ a well-mixed compartment assumption that ac-counts for transport by flow and diffusion using simplified mass-transfer coefficients.Other models fully account for spatial variations and transport to simulate smallvascular injuries under flow [24, 26] through the use of partial differential equations(PDEs). Despite the numerous mathematical models of flow-mediated intravascularblood clotting, there are relatively few models of hemostasis. To our knowledge, themodel we present here is the first mathematical model to incorporate key biophysicaland biochemical mechanisms of primary hemostasis in the framework of an extravas-cular injury. Our model is based on ODEs and is calibrated and validated usinganalogous experimental (Appendix A.2) and PDE (Appendix A.4) models [7, 40].One of the many contributions of this model is its computational efficiency, whichallows for systematic exploration of large sets of parameters that would be impossibleto resolve using the spatial PDE model or an experimental model. Parameters thatcontribute to bleeding include, but are not limited to, injury size, hemodynamics, ves-sel wall composition, and the anatomy of the adjacent perivascular and extravascularspaces. An experimental platform for studying hemostasis with a focus on collagen-TFinduced thrombus formation within a model vascular wall was previously presentedin [40]. This experimental model is an in vitro microfluidic flow assay (MFA) werefer to as the “bleeding chip”. Our mathematical model design is inspired by thespecific ‘H’-geometry that is characteristic of the bleeding chip and will be used inthe future to quantify variability of MFAs and identify components of the clottingsystem that underscore the known phenotypic variability in certain bleeding disor-ders [32]. The goal of this work is to create an efficient, mechanistic model that can This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. (A) (B) B
1. Blood Flow Through Injury2. Platelet Interactions3. ADP Dynamics4. Aggregate Growth1. Blood Flow Through Injury2. Platelet Interactions3. ADP Dynamics4. Aggregate Growth1. Blood Flow Through Injury2. Platelet Interactions3. ADP Dynamics4. Aggregate Growth1. Blood Flow2. Platelet Interactions3. ADP Dynamics4. Aggregate Growth
Injury Channel
Bottom Reaction ZoneTop Reaction Zone
Stokes Region A BloodWash
Fig. 1 . (A) Bleeding chip.
Blood and wash buffer are introduced into the two vertical chan-nels of the bleeding chip. A horizontal injury channel connects the two vertical channels, whereplatelet deposition occurs [40]. (B)
Mathematical model components.
1. Pressure driven flowwith a Stokes-Brinkman layered velocity calculation. 2. Platelets’ state (mobile/bound and unac-tivated/activated), mediated by activating chemical concentration [ADP], determines adhesion toinjury and cohesion to other platelets. (3) ADP is released from bound platelets and is subject totransport by both advection and diffusion. (4) Aggregate growth is influenced by changes in thicknessand porosity and feedbacks on blood flow due to increased frictional resistance.
2. Mathematical Model.
We model the situation where blood flows through ahole in the vessel wall into the extravascular space. In this scenario, a disruption in theendothelial layer occurs and exposes collagen to the blood. Platelets are transportedto the injury by flow and adhere to collagen. This in turn activates them, leading tothe release of platelet agonist ADP and the recruitment of additional platelets thatcohere and form an aggregate. The model incorporates the flow of blood throughan ‘H’-geometry to mimic the bleeding chip, as described above and as depicted inFigure 1A, and it explicitly models aggregate formation under local fluid dynamicconditions. Components of the model include descriptions of the blood flow throughthe bleeding chip geometry and the injury channel as well as the platelet interactions,soluble agonist ADP dynamics, and aggregate growth (see Figure 1B).Each component of the model involves a number of compartments which cor-respond to different spatial portions of the injury channel. Two compartments aredefined by the platelet aggregates on the bottom and top walls of the channel asshown in Figure 2A. Others are defined in terms of diffusive boundary layers forplatelets (Figure 2B) and the activating chemical ADP (Figure 2C). There are dif-fusive boundary layers associated with each of the aggregates. Lastly, there is acompartment representing the ‘gap’ between the boundary layers associated with the
This manuscript is for review purposes only.
K.G. LINK, ET AL. (A) (B) (C) A Top AggregateBottom Aggregate d agg d agg B Top Platelet Reaction ZoneBottom Platelet Reaction Zone d prz d prz C ADP in Bottom ADP in Top AggregateADP in Bottom AggregateADP in Top BLADP in Gap d abl d abl d gap Fig. 2 . Spatial compartments of model.
Each panel shows spatial compartments of theinjury channel. (A) Two aggregate (Brinkman) compartments with thickness d agg with a volumeof V agg and the Stokes region. (B) Two aggregate compartments and the platelet boundary layers(PBL), together defining platelet reaction zones (PRZ) with thicknesses of d prz and volumes V prz describe a two compartment platelet model. (C) Two aggregate compartments and the boundarylayers associated with ADP (ABL) with thicknesses d abl . The remaining compartment in the bulkof the Stokes region is denoted as the “gap” and has thickness d gap . These layers define a fivecompartment ADP model. bottom and top aggregates shown in Figure 2C. The two compartments of the plateletmodel and the five compartments describing the ADP model are defined precisely be-low. The model unknowns are number densities for three populations of platelets andthe concentration of ADP in the various compartments; the volume fraction of boundplatelets in and the thickness of each of the aggregate compartments, and the fluidvelocity profile across the channel. The velocity profile is the only fully spatially-dependent variable in the model; it is used in defining the diffusive boundary layersand in defining advective fluxes of fluid, platelets, and ADP into and out of the variouscompartments. The mathematical representation of flow in this geom-etry is nonstandard. Below, we describe how we model pressure-driven flow in theinjury channel and incorporate hinderance by the porous, growing aggregate. Theflow of blood through the bleeding chip is described in Subsection 2.1.1. Output fromthis calculation determines the inputs used in the Stokes-Brinkman calculation of theflow inside of the injury channel where aggregates form. The resulting velocity profileis used to calculate the increased resistance in the injury channel and therefore altersthe pressure drop P − P . Details of this calculation are found in Subsection 2.1.2. The bleeding chipcontains two vertical channels ( L × w × L z ) = (10 mm × µ m × µ m), con-sisting of a ‘blood’ channel and a ‘wash’ channel, connected by a horizontal ‘injury’channel ( L h × d h × L z ) = (150 µ m × µ m × µ m), which forms an ‘H’-shapedgeometry (Figure 1A). The flow rates and/or pressures at the inlets and outlets areprescribed and blood flows into the device through the blood channel. A fraction ofthe blood exits through the outlet of the blood channel. The remaining blood makesits way through the injury channel (from right to left) and exits through the outlet ofthe injury channel into the wash channel. Two hydraulic resistors in series describethe flow in each of the blood and wash channels and a connecting resistor in paralleldescribes the flow in the injury channel (Figure 3A). The resistances R , R , R , R are input parameters; the resistance R M is obtained from our calculation of the flowin the injury channel as explained in the next section. By solving the linear system ofequations for this hydraulic circuit (HC) as shown in Appendix A.1, we determine thepressure gradient G h = ( P − P ) /L h across the injury channel. The resulting pressuregradient G h is used to determine the flow velocity through the injury channel. This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. (A) (B) R R R y = 0 y = d y = d h d y = d h L h L z x yz TOPBOTTOM R top R bot R SR BA x yz R t R s R b L z d h L h B P B P CB P W P CW P P Q Q Q Q Q M R M R R R R R R R y = 0 y = d y = d h d y = d h L h L z x yz TOPBOTTOM R top R bot R SR BA x yz R t R s R b L z d h L h B P B P CB P W P CW P P Q Q Q Q Q M R M R R R R Fig. 3 . (A) Hydraulic circuit (HC) system describing flow through the bleeding chip. Quantitiesin blue are inputs while quantities in red are the desired output determined by solving system ofequations in Appendix
A.1 . (B) A schematic of the domain filled with Brinkman regions greenrepresenting the bottom R b and top aggregates R t and a Stokes region R s with a depth of L z , widthof d h , and a length of L z denoted in orange. Flow moves from right to left and the resulting velocityprofile is depicted in blue. The calculation of the fluid velocity profile ateach time proceeds as follows. We assume that we know the current thicknesses d b and d t of the aggregate compartments (see Figure 2A), and the volume fractions ofbound platelets θ B b and θ B t of the aggregates. Define the intervals R b = { ≤ y ≤ d b } , R s = { d b ≤ y ≤ d h − d t } , and R t = { d h − d t ≤ y ≤ d h } . Given that the Reynoldsnumber in the injury channel of the bleeding chip is R e ≈ . << x -direction for y ∈ R s is modeled by theStokes equations, and the flows for y ∈ R b and y ∈ R t are modeled by the Brinkmanequations, as we have done in previous studies [23, 24]. These are the Stokes equationsmodified by the inclusion of a Brinkman drag term among the forces acting on thefluid. The fluid dynamics equations are0 = − ∂p b ∂x + µ ∂ u b ∂y − µα b ( θ B b ) u b , for y ∈ R b , (2.1) 0 = − ∂p s ∂x + µ ∂ u s ∂y , for y ∈ R s , (2.2) 0 = − ∂p t ∂x + µ ∂ u t ∂y − µα t ( θ B t ) u t , for y ∈ R t . (2.3)where we assume the density of the fluid is ρ = 1 g/cm and the Brinkman coefficients α b and α t are functions of θ B b and θ B t , respectively. These differential equations aresupplemented with no-slip boundary conditions u b (0) = 0 and u t ( d h ) = 0 as wellas matching conditions at the edge of each aggregate compartment. The matchingconditions are that the velocity u and the shear stress µ dudy are continuous for all y .The pressure gradient ∂p∂x in the injury channel is assumed to be independent of y , soit is a constant G h determined from the circuit calculation found in Appendix A.1.To relate the volume fraction of bound platelets θ B b , θ B t to the permeability (1 /α b ,1 /α t ) of the Brinkman layers, the functions α b and α t are defined with the widely-used This manuscript is for review purposes only.
K.G. LINK, ET AL.
Kozeny-Carman relation [30]. α ( θ ) = C K θ (1 − θ ) , and C K is a fitted parameter whose estimation is discussed in Appendix C. Becauseof the linearity of the differential equations above (2.1)-(2.3), the velocity in eachcompartment is a linear combination of exponential functions of y , with a total of sixunknown coefficients. We determine these coefficients by solving the linear systemof equations corresponding to the two boundary conditions and the four matchingconditions (Appendix B). For later convenience, we define the velocity as u ( y ) =( u b ( y ) for y ∈ R b , u s ( y ) for y ∈ R s , u t ( y ) for y ∈ R t ).In practice, we first determine a preliminary velocity field by taking G h = 1.This allows us to calculate R M = L h · (cid:82) dh u ( y ) dy for use in the circuit analysis describedin Appendix A.1. Once the actual G h is known from the circuit analysis, the velocityis obtained by multiplying the preliminary velocity by G h . The velocity profile andthe upstream concentration of platelets dictate the transport of mobile platelets andADP to and from the site of injury. Platelets and ADP both move by advection anddiffusion. For ADP, diffusion is Brownian. For platelets, diffusion is used as a modelfor the effects on platelet motion of local flow disturbances generated by the complexmotions of the deformable red blood cells that make up approximately 40% of thebloods volume [42]. The diffusion of platelets and ADP to and from the walls/growingaggregates define the regions in which platelet activation, adhesion and cohesion canoccur. For clarity and brevity, all descriptions of the dynamics of platelets, solubleagonist ADP, and aggregate growth will be associated with the bottom wall of theinjury, although there are similar equations for dynamics near the top wall. Wetherefore will temporarily drop the subscript b that denotes ‘bottom’. We consider a three species model of platelet aggregation involv-ing the number densities (plt/ µ L) of unactivated and activated mobile platelets P m,u , P m,a in the platelet reaction zone (PRZ) and bound platelets P b,a in the growing ag-gregate region. The PRZ is the union of the aggregate region and the adjacent plateletboundary layer (PBL) (Figure 2B). Note the superscripts m and b are associated witha platelet’s state of mobility (mobile or bound) and the superscripts u and a indicateits state of activation (unactivated or activated). The first two platelet populationsare mobile and are subject to the effects of flow within and (sufficiently) near thegrowing aggregate. P m,u , enter upstream and leave downstream whereas P m,a onlyleave downstream because we have assumed that no mobile, activated platelets comeinto the injury channel from upstream. P m,u can adhere to the wall and thereforebecome bound and activated. Activated, mobile platelets can either adhere to thewall or cohere to P b,a . In addition to adhesion, P m,u are activated by exposure to thesoluble agonist ADP. Platelet adhesion and cohesion increase the volume of boundplatelets in the aggregate and consequently increase the thickness and volume fractionof bound platelets of each aggregate. They thus indirectly affect the delivery/removalof mobile platelets by flow. Similarly to what we did in[15, 22], we approximate the effects of diffusive motion in conjunction with advectivetransport by defining appropriate boundary layers, the regions from inside of whichplatelets can reach the growing aggregate before being carried away downstream. Weestimate the thickness of the boundary layer as follows (our estimates agree very well
This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. d agg and wedefine a platelet boundary layer as follows. For each x ∈ (0 , L h ), let h plt ( x ) be definedby equating the typical time it would take a platelet at location ( x, d agg + h plt ( x )) todiffuse to the edge of the aggregate with the time it takes that platelet to be carriedto the outlet of the injury by the flow: h plt ( x ) D p = L h − xu ( d agg + h plt ( x )) . Solving for h plt ( x ) defines a platelet boundary layer for platelets as shown in Fig-ure 4A). Because the platelet has a non-vanishing size, we modify our definition ofthe top edge of the platelet boundary layer to d pbl ( x ) = min( h plt ( x ) , d plt ), where d plt is the diameter ( ≈ µ m) of a platelet. We are interested in the advective flux ofplatelets into and out of this compartment, shown in Figure 4B. Transport into thiscompartment at the inlet of the injury channel depends on the inlet platelet concen-tration and the fluid velocity u ( y ) for d agg < y < d agg + d pbl (0). While the boundarylayer narrows between the inlet and the outlet of the injury channel, platelets canleave the boundary layer by advection for this same range of y values. Hence, letting d pbl = d pbl (0), we define the PBL compartment associated with the bottom aggregateto be those points in the injury channel with d agg < y < d agg + d pbl . For the bottomaggregate, the PRZ is defined by 0 ≤ y ≤ d prz = d agg + d pbl . A similar constructionis used to define the platelet boundary layer compartment adjacent to the aggregateon the top injury channel wall. (A) (B) d agg h plt ( x ) A FLOW ( x, d agg + h plt ( x )) d agg B FLOW d prz = d agg + d pbl Fig. 4 . Boundary layers associated with the transport of platelets.
Blood flows fromright to left, entering the injury channel at x = 0 and exiting downstream at x = L h . (A) Consider aplatelet at the position ( x, d agg + h plt ( x )) . The distance it must travel to reach the edge of the growingaggregate is h plt ( x ) (cyan) and the thickness of the aggregate is denoted as d agg . (B) The PBLthickness at the entrance of the injury channel h plt (0) (cyan, dashed) defines the reaction zone forplatelet interactions in that d prz = d agg + d pbl , where d pbl = min( h plt (0) , d se ) and d se = d plt = 2 µ m. The advective transport of platelets into and out of each of the platelet compart-ments is hindered to a degree that depends on the fraction of each compartment’svolume that is occupied by platelets. Let θ M = v plt ( P m,u + P m,a ) and θ B = v plt P b,a ,where v plt is the volume of an individual platelet. The quantities θ M and θ B arethe volume fractions occupied by mobile and bound platelets, respectively. Note thatbound platelets are found only in the aggregate compartments while mobile plateletscan be found throughout the injury channel, but we track separately the numberdensities of mobile platelets in the two PRZs. To describe the transport of mobile, This manuscript is for review purposes only.
K.G. LINK, ET AL. unactivated platelets into the injury channel, we define the quantity J in prz = (cid:16) (cid:90) d prz d agg u ( y ) P up ( y ) dy (cid:17) W ( θ M ) + (cid:16) (cid:90) d agg u ( y ) P up ( y ) dy (cid:17) W ( θ M + θ B ) , where P up ( y ) is the number density of mobile, unactivated platelets at the inlet to theinjury channel. In this work, we assume a uniform and constant upstream distribution P up ( y ) = P up, ∗ , but using non-uniform distributions is straightforward (Appendix D).Note that the use of P up, ∗ in the model yields symmetric top and bottom aggregates.The model framework allows for asymmetric aggregate formation if P up ( y ) is notsymmetric and simulations found in Appendix D explore the resulting aggregates.The function W is intended to limit platelet entry and movement through the regionin which the volume fraction of platelets is already high. We used the specific function W ( θ ) = tanh( π (1 − ( θ/θ max ))) , which monotonically decreases from W (0) = 1 to W (1) = 0 [24]. Therefore, thereis no hinderance when the volume fraction is 0 and total hinderance as the volumefraction approaches 1. In choosing this function, we assume that the ability for mobileplatelets to move into a region is only gradually impaired until θ nears θ max = 0 . L z J in prz .We define the following quantities to describe the transport of unactivated andactivated mobile platelets out of the injury channel. The rates at which these plateletsleave the PRZ are L z J out,u prz and L z J out,a prz , respectively, where J out,u prz = (cid:34)(cid:16) (cid:90) d prz d agg u ( y ) dy (cid:17) W ( θ M ) + (cid:16) (cid:90) d agg u ( y ) dy (cid:17) W ( θ M + θ B ) (cid:35) P m,u ,J out,a prz = (cid:34)(cid:16) (cid:90) d prz d agg u ( y ) dy (cid:17) W ( θ M ) + (cid:16) (cid:90) d agg u ( y ) dy (cid:17) W ( θ M + θ B ) (cid:35) P m,a . The flow of blood in the injury channel is coupled to the transport of platelets intoand out of the injury channel as we describe below.
Platelets contribute to aggregateformation through two binding processes: adhesion and cohesion. Both unactivatedand activated mobile platelets can adhere directly to the bottom or top wall of theinjury channel. Mobile, activated platelets can cohere to already bound, activatedplatelets. Platelets are activated through direct contact with the wall or by exposureto the soluble agonist ADP. The activation state of the platelet as well as the availablespace both on the wall and in the developing aggregate determine which bindingprocesses can occur.Since adhesion requires that a platelet directly contacts the wall, we limit adhesionto those platelets within a specified distance d se from the walls. In the model, aplatelet’s location relative to the walls is only known with respect to whether theplatelet is inside one of the PRZs, and so we enforce this limit approximately in theadhesion rate function k adh ( d agg , θ B ) = k adh min (cid:18) − d agg d se , (cid:19)(cid:18) θ max − θ B (cid:19) , This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. k adh is a first order rate constant, the second term represents the accessibilityof the wall, and the third term represents the porosity of the aggregate. With thisfunction, adhesion is further limited as the volume fraction of platelets in the aggregateapproaches its maximum value θ max , because platelets already bound to the walloccupy a portion of the wall space available for adhesion. Both unactivated andactivated mobile platelets can adhere to the wall and hence the rate of adhesion isproportional to the adhesion rate function and the total number density of mobileplatelets.For mobile platelets to cohere to the aggregate, they must be activated and suffi-ciently close to activated bound platelets. As previously mentioned, activated plateletshave activated integrin receptors on their surfaces that allow for α IIb β -fibrinogen- α IIb β and α IIb β -vWF- α IIb β bond formation. These bonds are strong, and long-lived with very small off rates. Therefore, we assume that bond breakage is negligibleand do not allow for platelet detachment in this model. The rate of platelet-plateletcohesion, k coh P b,a P m,a , depends on the second order rate constant k coh and the num-ber densities of bound activated platelets, P b,a , and mobile, activated platelets P m,a ,respectively.In addition to changing the number of bound platelets in the aggregate and theflow through the injury channel, the processes of adhesion and cohesion change thetotal volume of the aggregate and of the adjacent boundary layers (and therefore ofthe PRZs) by increasing their thicknesses. For a given number of platelets in thePRZ, an increase in that compartment’s volume results in a decrease in the numberdensity of platelets in it – an effect we refer to as dilution. For example, the termdescribing this effect on mobile unactivated platelets has the form (cid:32) ddt ( d prz ) d prz (cid:33) P m,u , and appears naturally in calculating the rate at which the number of these plateletsin the platelet reaction zone changes. Here we account for the processesdescribed in Subsection 2.2.1 and Subsection 2.2.2 to derive evolution equations forthe platelet number densities. Let V prz = L h L z d prz denote the volume of the plateletreaction zone. Then the rate of change of the number of mobile, activated plateletsin that reaction zone is d ( V prz P m,u ) dt (cid:124) (cid:123)(cid:122) (cid:125) rate of change of = L z ( J in prz − J out,u prz ) (cid:124) (cid:123)(cid:122) (cid:125) flow in/out of reaction zone − V prz k act P m,u (cid:124) (cid:123)(cid:122) (cid:125) activation − V prz k adh P m,u (cid:124) (cid:123)(cid:122) (cid:125) adhesion Applying the product rule on the left-hand side and rearranging terms, we determinethat dP m,u dt = ( J in prz − J out,u prz ) L h d prz (cid:124) (cid:123)(cid:122) (cid:125) flow in/out of reaction zone − k act P m,u (cid:124) (cid:123)(cid:122) (cid:125) activation − k adh P m,u (cid:124) (cid:123)(cid:122) (cid:125) adhesion − (cid:32) ddt ( d prz ) d prz (cid:33) P m,u (cid:124) (cid:123)(cid:122) (cid:125) dilution . (2.4)The first and third terms describe transport and adhesion. The second term describesactivation of mobile unactivated platelets by exposure to ADP. Based on similar This manuscript is for review purposes only. K.G. LINK, ET AL. considerations, we find evolution equations for the number densities of the otherplatelet species P m,a and P b,a . dP m,a dt = − J out,a prz L h d prz (cid:124) (cid:123)(cid:122) (cid:125) flow out of reaction zone + k act P m,u (cid:124) (cid:123)(cid:122) (cid:125) activation − k adh P m,a (cid:124) (cid:123)(cid:122) (cid:125) adhesion (2.5) − k coh P b,a P m,a (cid:124) (cid:123)(cid:122) (cid:125) cohesion − (cid:32) ddt ( d prz ) d prz (cid:33) P m,a (cid:124) (cid:123)(cid:122) (cid:125) dilution ,dP b,a dt = d prz d agg (cid:16) k adh ( P m,u + P m,a ) + k coh P b,a P m,a (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) adhesion and cohesion − (cid:32) ddt ( d agg ) d agg (cid:33) P b,a (cid:124) (cid:123)(cid:122) (cid:125) dilution . (2.6) An unactivated platelet stores quantities of ADPin its dense granules and releases this ADP into the surrounding fluid when the plateletis activated. Similar to the treatment in [24], we assume this release occurs duringthe time interval 1-5s after activation. Once released, that ADP moves by advectionand diffusion within and out of the aggregate where it can activate mobile platelets.Describing these processes requires defining two boundary layers: one within the ag-gregate h adp agg ( x ) and one in the Stokes region h adp BL ( x ). These boundary layers are usedto define the ADP compartments shown in Figure 2C and the ADP concentrationstracked are those inside the aggregate [ ADP ] agg , in the Stokes region boundary layer[ ADP ] BL , and in the remaining Stokes region which we denote as [ ADP ] gap . (A) (B) d agg A FLOW ( x, d agg h adp agg ( x )) h adp agg ( x ) B FLOW ¯ h adp BL h adp BL (0)¯ h adp agg Fig. 5 . Boundary layers associated with the transport of soluble agonist ADP. (A)Consider an ADP molecule at the position ( x, d agg − h adp agg ( x )) . The distance it must travel to reachthe edge of the growing aggregate is h adp agg ( x ) . (B) The boundary layer thickness at the entrance of theinjury channel in the Stokes region is h adpBL (0) (magenta, dashed) and defines the reaction zone forADP interactions in that region (ARZ). Both the average thickness in the Stokes region ¯ h adpBL (blue,dashed) and in the growing aggregate ¯ h adp agg (purple, dashed) are used to define diffusive transportfrom the aggregate into the Stokes region. The first boundary layer is related toADP movement in the Stokes region after it has left the aggregate where it is released.The calculation is like that for platelets in the Stokes region as described above inSubsection 2.2.1. This boundary layer is thinnest at the upstream end and widestat the downstream end. Its thickness is greater than that, h plt ( x ), of the plateletboundary layer. This is because (i) ADP has a larger diffusivity than platelets and (ii)the choice of velocity used in the boundary layer calculation. For the platelets, it is avelocity in the Stokes region u ( d agg + h plt ( x )) and for the ADP it is the (lower) velocity This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. u ( d agg ). Therefore, h adp BL ( x ) > h plt ( x ) , ∀ x ∈ [0 , L h ].As shown in Figure 5B, the maximum thickness of this ADP boundary layer in theStokes region is denoted as h adp BL (0) and the average boundary layer thickness as ¯ h adp BL .The second boundary layer describes ADP movement within the aggregate to-wards the Stokes region and requires the new calculation given in this section. Wedo not have a separate ADP concentration in this boundary layer; its thickness isused only in calculating the ADP concentration gradient. Consider the starting point( x, d agg − h adp agg ( x )) inside the aggregate and equate the time it takes an ADP moleculethere to diffuse into the Stokes region with the time it would take it to be washedpast the injured region by the flow:( d agg − h adp agg ( x )) D a = L h − xu ( d agg ) . We solve for the thickness h adp agg ( x ) and define ¯ h adp agg to be the average aggregate layerthickness shown in Figure 5B. The ADP boundary layer thickness d abl = h adp BL (0) andthe thickness of the growing aggregate contribute to the total thickness of the regionsin which we track the ADP concentrations: [ ADP ] agg , [ ADP ] BL , and [ ADP ] gap . The source of ADP in this model cor-responds to the release of ADP from activated, bound platelets. The rate of releaseis σ release ( t ) = (cid:90) ∞ ˆ AR ( τ ) ddt ( P b,a )( t − τ ) dτ, where ˆ A (Appendix C) is the total quantity of ADP released by an activated platelet,ˆ AR ( τ ) is the rate of release at an elapsed time since activation τ and (cid:82) ∞ R ( τ ) dτ = 1.The function R ( τ ) is defined as having value zero up to one second after activation,a positive bell-shaped function for the time interval 1 < τ < ddt ( P b,a )( t − τ ) dτ is the number of platelets newly activatedand bound in the time interval [ t − τ, t − τ + dτ ]. We account only for ADP releasefrom bound activated platelets because activated platelet that do not become boundare carried downstream before they begin to release their ADP.To facilitate our describing the advective transport of ADP out of the aggregate,the boundary layer, and the gap, we define the expressions J adp agg = (cid:16) (cid:90) d agg u ( y ) dy (cid:17) [ ADP ] agg ,J adp BL = (cid:16) (cid:90) d agg + d abl d agg u ( y ) dy (cid:17) [ ADP ] BL ,J adp gap = (cid:16) (cid:90) d agg + d abl + d gap d agg + d abl u ( y ) dy (cid:17) [ ADP ] gap , respectively. Then, the rates of advective transport of ADP out of the ADP com-partments in the injury channel are L z J adp agg , L z J adp BL , and L z J adp gap . Early on, themagnitude of the velocity and the thicknesses of the aggregate and boundary layerslead to the ordering J adp gap > J adp BL > J adp agg . As the aggregate grows, the boundarylayers become larger, decreasing the thickness of the gap d gap and resulting in a newordering J adp BL > J adp gap . If the flow is slowed sufficiently, d gap = 0 and J adp gap = 0. This manuscript is for review purposes only. K.G. LINK, ET AL. (A) (B) J ⇤ di↵ J di↵ A FLOW J di↵ J ⇤ di↵ B FLOW ` di↵ = ¯ h adp BL + ¯ h adp agg ` ⇤ di↵ = ¯ h adp BL + d gap / Fig. 6 . Compartment model describing diffusive transport of ADP. (A) Diffusivefluxes of ADP from the aggregates (green) to the associated boundary layers (blue) are denoted as J diff . ADP also diffuses from the boundary layers associated with the bottom and top aggregatesto the bulk of the Stokes region (pink) with flux J ∗ diff . The fluxes are subject to changing layerthickness that are determined by the flow velocity. (B) The dashed lines correspond to boundarylayers associated with the bottom aggregate. The specified lengths (cid:96) diff and (cid:96) ∗ diff are used to definethe gradients of ADP concentration between the aggregates and boundary layers as well as betweenboundary layers and bulk Stokes region, respectively. Diffusive transport of ADP from the ag-gregate into the boundary layer and into the bulk Stokes region also changes the ADPconcentrations. The diffusive transport of ADP (Figure 6A) between the aggregate(green) and the boundary layer region (pink) and that between boundary layer and thebulk Stokes region (grey) are determined using Fick’s law [11]. Let V abl = L h L z d abl be the volume of the ADP boundary layer and (cid:96) diff be the length that determines thegradient of ADP from the aggregate to the boundary layer. As shown in Figure 6B,this length is the distance between the edge of the boundary layer in the aggregate(purple, dashed) and the edge of the average boundary layer in the Stokes region(dashed, blue). The diffusive flux of ADP from the aggregate into the boundary layeris J diff = D a ([ ADP ] agg − [ ADP ] BL ) c (cid:96) diff , and the rate of diffusion is L h L z J diff , where L z L h is the area of the interface betweenthe zones, D a is the diffusion coefficient, ([ ADP ] agg − [ ADP ] BL ) is the difference inADP concentration between the aggregate and the boundary layer, and (cid:96) diff is thelength used to define the gradient. Similarly, the diffusive flux of ADP from theboundary layer into the ‘gap’ is J ∗ diff = D a ([ ADP ] BL − [ ADP ] gap ) c (cid:96) ∗ diff , where the rate of transport flux is L h L z J ∗ diff , ([ ADP ] BL − [ ADP ] gap ) is the differencein ADP concentration and (cid:96) ∗ diff is half the distance between the edge of the boundarylayer and the center of the gap. The fitted parameters c and c are discussed inAppendix C. Accounting for release, transport dueto advection, and diffusive transport between the aggregates, boundary layer regions,
This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. d [ ADP ] agg dt = σ rel ( t ) (cid:124) (cid:123)(cid:122) (cid:125) release − J diff d agg (cid:124) (cid:123)(cid:122) (cid:125) diffusioninto BL − J adp agg L h d agg (cid:124) (cid:123)(cid:122) (cid:125) flow out of agg − (cid:18) d ( d agg ) dt d agg (cid:19) [ ADP ] agg (cid:124) (cid:123)(cid:122) (cid:125) dilution , (2.7)and d [ ADP ] BL dt = J diff d abl (cid:124)(cid:123)(cid:122)(cid:125) diffusioninto BL − J ∗ diff d abl (cid:124)(cid:123)(cid:122)(cid:125) diffusioninto gap − J adp BL L h d abl (cid:124) (cid:123)(cid:122) (cid:125) flow out of BL − (cid:18) d ( d abl ) dt d abl (cid:19) [ ADP ] BL (cid:124) (cid:123)(cid:122) (cid:125) dilution , (2.8)The concentration of ADP in the gap changes because of diffusive transport from both the bottom and top boundary layers: d [ ADP ] gap dt = J ∗ diff,b d gap (cid:124) (cid:123)(cid:122) (cid:125) diffusion intogap from bot-tom + J ∗ diff,t d gap (cid:124) (cid:123)(cid:122) (cid:125) diffusion into gapfrom top − J adp gap L h d gap (cid:124) (cid:123)(cid:122) (cid:125) flow out of gap − (cid:18) d ( d gap ) dt d gap (cid:19) [ ADP ] gap (cid:124) (cid:123)(cid:122) (cid:125) dilution , (2.9)where J ∗ diff,b and J ∗ diff,t are the respective diffusive fluxes of ADP into the gap fromthe bottom and top boundary layers. The formation of the aggregate layer is describedby its evolving thickness, d agg , and volume fraction of bound platelets, θ B . To de-termine evolution equations for these quantities, one must consider the adhesion andcohesion events that contribute to the change in volume occupied by bound platelets.Let V agg = L h L z d agg . The rate of change of the volume occupied by bound plateletsis ddt ( V agg θ B ) = V prz k adh ( P m,u + P m,a ) v plt (cid:124) (cid:123)(cid:122) (cid:125) adhesion of mobile plts. + V prz ( k coh P b,a P m,a ) v plt (cid:124) (cid:123)(cid:122) (cid:125) cohesion of mobile, act. plts . Dividing by d agg θ B and using the product rule, and partitioning the adhesive termand cohesion term we find the following equations. dθ B dt = γ (cid:18) d prz d agg (cid:19) k adh (cid:16) P m,u + P m,a (cid:17) v plt + β (cid:18) d prz d agg (cid:19) k coh (cid:16) P b,a P m,a (cid:17) v plt (2.10) ddt ( d agg ) = (1 − γ )( d prz ) k adh (cid:18) P m,u + P m,a P b,a (cid:19) + (1 − β )( d prz ) k coh (cid:18) P b,a P m,a P b,a (cid:19) , (2.11)where β , γ are assumed to have the functional forms γ = (cid:40) , if θ B = θ max , if d agg ≥ d plt β = (1 + tanh( k b ( s b − θ B ))) , if d gap > , θ B < θ max , if θ B = θ max , if d gap = 0 . This manuscript is for review purposes only. K.G. LINK, ET AL.
The partition functions β and γ are phenomenological in nature and will be par-tially validated by PDE model (Appendix A.4). Details are found in Appendix C.Equations (2.4)-(2.11) describe platelet dynamics, ADP dynamics, and aggregate for-mation associated with the bottom wall of the injury zone. There is a similar set ofequations governing the concentrations of platelets species, ADP, aggregate thickness,and aggregate volume fraction associated with the top aggregate. As noted above, once an unactivated platelet entersthe PRZ it can adhere to the wall and therefore become bound and activated, it can beactivated by ADP and possibly cohere to already bound platelets, or it can be carriedout of the injury channel by the flow. While platelet activation is a complex processtriggered by a range of stimuli and made up of an ensemble of responses [27], we focuson activation by direct adhesion to the wall and by exposure to ADP and on theresulting activation of integrin α IIb β receptors and release of ADP. We discussed theactivation by direct adhesion above. Activation triggered when ADP in the plasmabinds to a platelet’s P2Y and P2Y receptors, provides a means of activation thatdoes not require contact with the SE matrix. Because ADP stimulates this activationand additional ADP is released as a consequence, there is a positive feedback in thesystem involving the concentration of ADP and the transition of mobile, unactivatedplatelets to mobile, activated platelets. To model ADP-dependent activation, weassume that the activation trigger is a weighted average of the ADP concentrationsin the aggregate and in the boundary layer. We define the rate of ADP-inducedactivation as k act = k act A ([ ADP ] agg , [ ADP ] BL ) , where A ([ADP] agg , [ADP] BL ) = (cid:18) d agg d arz (cid:19) [ADP] agg [ADP] ∗ + [ADP] agg + (cid:18) d pbl d arz (cid:19) [ADP] BL [ADP] ∗ + [ADP] BL . Here, d arz = d agg + d abl is the thickness of the ADP reaction zone (ARZ) and [ ADP ] ∗ =1 µ M. For some of our simulations, we instead considered ADP-independent plateletactivation in the plasma by setting A ([ADP] agg , [ADP] BL ) ≡
1, in which case onlythe value of the rate constant k act matters. Below we explore the effects of ADP-independent versus ADP-dependent activation on aggregate growth and occlusion ofthe injury channel. Simulations of the ODE model consist of solvingthe HC system and the Brinkman-Stokes system to determine G h and R M , respec-tively. The differential equations describing the evolution of platelet species, ADP, andaggregate thickness and density are solved using MATLAB function ode15s , whichemploys a variable-step, variable-order (VSVO) solver based on the numerical differ-entiation formulas (NDFs) of orders 1 to 5 for stiff systems [41]. During each timestep of the computation, we perform the following series of updates for the unknowns.1. Assume G h = 1 and solve the linear system in Appendix B to determine u ( y )and R M .2. Use R M as an input parameter for the HC system in Appendix A and solvefor G h .3. Scale the flow velocity using the pressure gradient update; ¯ u ( y ) := G h u ( y ).4. Use u ( y ) and diffusivity constants ( D p , D a ) to determine boundary layers( h plt ( x ), h adp BL ( x ), h adp agg ( x )) for both aggregates. This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. d prz ), compartments in the Stokes regiontracking ADP ( d abl , d gap ), and the growing aggregate ( d agg ) associated withthe bottom and top walls.6. Count the platelets activated within the previous time step and update theADP release function σ rel for each aggregate.7. Update mobile platelets and ADP concentrations associated with the topand bottom aggregates and reaction zones to account for advection using thevelocity profile u ( y ).8. Update ADP concentrations to account for diffusion using ADP boundarylayer thicknesses (¯ h adp BL , ¯ h adp agg ) of the top and bottom ARZs and aggregates.9. Update the thicknesses of the aggregates d agg , the volume fractions of boundplatelets θ B , and hinderance function W ( θ B + θ M ) to account for plateletadhesion and cohesion events.10. Repeated until desired time.The simulations described below were carried out for 20 minutes and all parametersvalues are found in Table 1-Table 7 with specific details in Appendix A - Appendix C.Each simulation requires ≈
50 seconds of computational time. This computationalcost can be substantially reduced by using a compiled language.
3. Results.
The presented mathematical model was calibrated and validatedthrough comparisons of model output with that of analogous MFA (Appendix A.2)and PDE (Appendix A.4) models. Preliminary studies with the PDE model revealedthat the pressure drop ∆ P = P − P across the injury channel was not constant withrespect to the y-direction. Therefore to achieve the most accurate model comparisons,it was necessary to compute an approximation of the flow rate Q P DE through theinjury channel using the PDE model for all the following simulations. This changesthe second step of our numerical scheme. Specifically, we do not solve the HC systemfor the pressure gradient G h . Instead, we assume G h = Q P DE ( d b , d t , θ Bb , θ Bt ) R M L eff h , where d b , d t and θ Bb , θ Bt are the thicknesses and densities of the bottom and top aggre-gates, respectively. The effective length of the channel, L eff h = 165 µ m as described inAppendix A.1. Details regarding Q P DE are found in Figure 7A-B and Appendix A.3.Note that the ODE model is not dependent on the use of the PDE flow map and theuse of the ODE model beyond the scope of this work is independent of PDE modelresults.Using this new calculation for the pressure gradient, we validated the fluid calcu-lation associated with the ODE model in Figure 7C-D. Assuming ADP-independentplatelet activation ( A ([ ADP ] agg , [ ADP ] BL ) = 1), we used MFA occlusion times to de-termine physiologically relevant values for the kinetic rate constant k adp (Figure 8A-B). Fixing a value for k adp , we compared occlusion metrics Figure 8C-D from the ODEand PDE models. Next, we assessed and characterized occlusion of the injury channelas a function of the flow shear rate and k adp in Figure 9A-B and then conducted sim-ilar experiments to test the sensitivity of model output to ADP-dependent plateletactivation as shown in Figure 10A-D. Evidence of the effects of dilution of ADP dueto advective and diffusive transport on aggregate formation is found in Figure 11A-H. For this study, we simulated 2Dflow through the injury channel with fixed pre-formed platelet aggregates. For each
This manuscript is for review purposes only. K.G. LINK, ET AL. (A) (B) A Aggregate Density A gg r ega t e T h i ck ne ss B Relative Difference
Aggregate Density (C) (D)
C FLOW
D FLOW
Fig. 7 . Model calibration and validation. (A) 3D flow rates Q PDE in prebound plateletsimulations with the PDE model using spatial discretization ∆ x = 2 µ m and (B) the relative differ-ence between the flow rates generated with ( ∆ x = 2 µ m) and ( ∆ x = 1 µ m). (C) Velocity profiles fromODE and PDE models for an asymmetric prebound platelet configuration ( d t = 6 µ m, θ Bt = 0 . , d b = 8 µ m, θ Bb = 0 . ). (D) Similar prebound platelet configuration ( d t = 12 µ m, θ Bt = 0 . , d b = 4 µ m, θ Bb = 0 . ) experiment tracks the amount of ADP in the aggregate for the PDE andODE models. Note θ max = 0 . . simulation, a different symmetric configuration of preformed aggregates was specifiedby choosing the thicknesses d b = d t ∈ { , , , , µm } and bound platelet densities θ B b = θ B t ∈ [0 . , . , · · · . Q P DE for each configurationas detailed in Appendix A.3. The ‘flow map’ Q P DE ( d b , θ B b , d t , θ B t ) obtained when thePDE calculation used a spatial discretization with step size ∆ x = 2 µ m is shown inFigure 7A, and Figure 7B presents evidence that the flow rate calculation is insensitiveto further refinement in ∆ x . The maximum relative difference in flow rates for the twostepsizes, which occured for an aggregate density of ≈ . ≈ . µm , was approximately 8(10) − .Figure 7C shows velocity profiles from the ODE and PDE models for asymmetricaggregates. The ODE profile shows excellent agreement in both Brinkman layers and This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. (A) (B)(C) (D) Fig. 8 . Model calibration and validation. (A) Occlusion times generated by the MFA bleed-ing chip and (B) ODE model times to Q min = 0 . µ L/min as a function of the ADP-independentplatelet activation rate k act . Comparison of ODE (orange) and PDE (blue) model outputs witha fixed shear rate ( ≈ s − ) and ADP-independent platelet activation rate ( k act = 3 . s − ).(C) Total number of bound platelets in the aggregate and (D) rate of blood flow through the injurychannel. the Stokes region with the PDE profile obtained half-way through the injury channel.The flow through the top aggregate, which is thin ( d t = 4 µ m) but dense ( θ B t = 0 . d b = 8 µ m) yet loose ( θ B b = 0 .
06) bottom aggregate,because resistance to flow increases as aggregate density increases.In Figure 7D, we show results from another asymmetric pre-formed aggregateconfiguration, assuming that all of the platelets in the aggregates became activated attime t = 0. Shown are the amounts of ADP in each of the aggregates due to releaseof ADP by these platelets in both the ODE and PDE models. Less ADP is washeddownstream from the thin and dense bottom aggregate because the flow in it is slower.Consequently, the maximum concentration of ADP in the bottom aggregate is higherthan that in the thick and loose top aggregate (not shown). It is important to notethat the concentration of ADP in each aggregate is subject to diffusive transport.The ODE model accounts for lateral diffusion (perpendicular to the channel walls)but not longitudinal diffusion. The differences in the models’ outputs are largely dueto longitudinal diffusion of ADP out of the injury channel in the PDE model (resultsnot shown).One of our main objectives in this work was to develop a mathematical model This manuscript is for review purposes only. K.G. LINK, ET AL. that can generate aggregates that grow to fill the injury channel. This objective wasmotivated by experiments conducted in the microfluidic bleeding chip model (MFA).Figure 8A shows three metrics of occlusion times from experiments using the MFA.The median times for all three metrics of occlusion are approximately 3.25 minutes,with a minimum occlusion time of 1.15 minutes and a maximum time of 6.25 minutes.The metric of occlusion best suited for comparison with the ODE model output isthe ‘time to Q min = 0 . µ L/min measurement’. Using an ADP-independent plateletactivation regime ( A ([ADP] agg , [ADP] BL ) ≡
1) in the ODE model simulations yieldsthe ‘time to Q min ’ results obtained for a range of platelet activation rates as shown inFigure 8B. For k act > . − , the time to Q min is within the bounds defined by theMFA results. In Figure 8C-D, we compare the growth of aggregates using the ODEmodel (orange) with output from the PDE model (blue). Metrics are the timecourseof the total number of bound platelets and of the flow rate through the injury channel.For a fixed shear rate ≈ − and a fixed platelet activation rate of 3 . − , modelresults fall within the range of the experimental values. Note that in Figure 8B, thetime to Q min with the chosen activation rate is within the range (black, dashed lines)seen with the MFA. Additionally, both computational models yielded aggregates thatreduced the flow rate of the injury channel in an amount of time similar to that seenin the MFA results. To betterunderstand how shear rate and ADP-independent platelet activation affect aggregateformation, we performed simulations in which we varied both the shear rate andactivation rate constant k act . Aggregate thickness d agg and aggregate density θ B after (A) (B) Aggregate Thickness (microns)
Activation Rate (1/s) S hea r R a t e ( / s ) Aggregate Density
Activation Rate (1/s) S hea r R a t e ( / s ) Fig. 9 . Occlusion in the Injury Channel with ADP-independent activation.
ODEmodel aggregate (A) thicknesses ( µ m) and (B) densities due to variations in ADP-independentplatelet activation rates k adpact and shear rate γ , after 10 minutes of simulation time.
10 minutes of simulation time are shown in Figure 9A-B. For shear rates 200 s − <γ < − , occlusive aggregates formed if k act > . − while significantly thinnerbut high density aggregates form for k act below approximately 1.5 s − . For shear ratesless than 100 s − , aggregate thickness increases as k act is increased, but in contrast towhat is seen for high shear rates, the aggregate density decreased as the activation rateconstant was increased from 1 . . − . Hence the model predicts that occlusiveaggregates form under a wide range of shear rates provided platelet activation rate islarge enough, but the nature of the aggregates is sensitive to activation rate especiallyat low shear rates. This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. In this section we turn to the model’sbehavior when the platelet activation rate is dependent on the ADP concentration k pltact = k act A ([ADP] agg , [ADP] BL ). If the ADP concentration is substantially belowthe transition concentration [ADP] ∗ = 1 µ M in this function, results might be verydifferent from those just reported with a constant rate of platelet activation. (A) (B)
Activation Rate (1/s) S hea r R a t e ( / s ) Peclet Number [ADP] in Reaction Zone
Activation Rate (1/s) S hea r R a t e ( / s ) (C) (D) Aggregate Thickness (microns)
Activation Rate (1/s) S hea r R a t e ( / s ) Aggregate Density
Activation Rate (1/s) S hea r R a t e ( / s ) Fig. 10 . Aggregate Heterogeneity under variations in flow conditions and ADP-dependent platelet activation. (A) Peclet number describing the advective versus diffusive trans-port of soluble agonist ADP. (B) Concentration of ADP in the aggregate and associated boundarylayer. (C) Aggregate thickness and (D) aggregate density as functions of shear rate and plateletactivation rate constant. These results identified two flow regimes for a fixed platelet activation rateconstant k act = 3 . s − ; high shear (cyan square) and low shear (orange square). The metrics aretaken after 10 minutes of simulation time. We explored how different forms of transport within the aggregates affects theirgrowth by calculating the ADP-related Peclet number after 10 minutes of simulationtime. This non-dimensional number is defined as the ratio of the rate of advection byflow and the rate of diffusion of ADP. For our purposesPe = d abl u ( d apz ) D a , where u ( d apz ) is the velocity at the edge of the ADP boundary layer in the Stokesregion ( d apz = d agg + d abl ). Figure 10A shows that advective transport dominatesfor shear rates greater than 500 s − for activation rate constants k act ranging from0 . − . − , but that, while still large, it is much lower for shear rates less than This manuscript is for review purposes only. K.G. LINK, ET AL.
Time (min.) [ A D P ] ( m u M ) [ADP] at High Shear Aggregate (agg)Boundary Layer (BL)Stokes Region (gap)
Time (min.) [ A D P ] ( M ) [ADP] at Low Shear Aggregate (agg)Boundary Layer (BL)
Distance from Bottom Wall (microns) V e l o c i t y ( c m / s ) Velocity Profiles at High Shear
Distance from Bottom Wall (microns) V e l o c i t y ( c m / s ) Velocity Profiles at Low Shear
Time (min.) T h i ck ne ss ( m i c r on s ) Layer Thicknesses at High Shear
Aggregate (agg)Boundary Layer (BL)Stokes Region (gap)
Time (min.) T h i ck ne ss ( m i c r on s ) Layer Thicknesses at Low Shear
Aggregate (agg)Boundary Layer (BL)Stokes Region (gap)
Time (min.) R a t e ( p l t/ s ) Rates of Aggregate Growth - HS
Adhesion of P m,u
Adhesion of P m,a
Cohesion of P m,a
Time (min.) R a t e ( p l t/ s ) Rates of Aggregate Growth - LS
Adhesion of P m,u
Adhesion of P m,a
Cohesion of P m,a
Fig. 11 . Effect of ADP Dilution on Aggregate Formation.
ADP concentrations in theaggregate, boundary layer, and gap under (A) high and (B) low shear conditions. Note there is nogap in the low shear conditions due to the thicknesses of the aggregate and boundary layer. Thegrey-scale squares denote the times at which the velocity profiles were determined. Velocity profilesat 0.3, 5, 8.75, and 15 minutes under (C) high and (D) low shear conditions. Thicknesses of theaggregate, boundary layer, and gap under (E) high and (F) low shear conditions. The rates ofadhesion and cohesion under high (G) and low (H) shear conditions.
This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. − . We can see that for most of the shear rate range, the concentration of ADPin the PRZ after 10 minutes is much less that the activation transition concentration[ADP] ∗ = 1 µ M and that the shear rate must be reduced to below ≈
100 s − for theADP concentration to reach 0 . µ M. As shown in Figure 10B, increasing the valueof k act results in higher ADP concentrations. Plots of aggregate thicknesses anddensities in Figure 10C-D identify two flow regimes that for a fixed activation rateconstant yield occlusive aggregates under low shear conditions and thin non-occlusive,but dense, aggregates under high shear conditions.To better understand the differences in aggregates that form under high shear(cyan squares in Figure 10) and low shear (orange squares in Figure 10) conditions,we consider the model output in Figure 11A-H. The initial rate of adhesion is greaterin high shear case (Figure 11G,H) because of the faster replenishment of platelets bythe faster flow. Hence, the initial rate of ADP release is also greater in this case.Despite this, the ADP concentrations are more than 6-fold lower (Figure 11A,B)because the faster flow carries ADP away more rapidly. In the low shear case, thetwo ADP boundary layers (ABL) are so thick (Figure 11E,F) that there is no gapregion in that case. For both shear conditions, the velocity through the injury channeldecreases as the aggregates form. The reduction in peak velocity is approximately30% after 15 minutes in the high shear case (Figure 11C), while the relative reductionis much greater, at approximately 80%, under low shear (Figure 11D). This largerreduction in velocity results in lower rate of advective transport and consequentiallyincreases the boundary layer thickness (orange, dashed) for ADP (Figure 11F). Asthe aggregate grows, this layer becomes smaller until the injury channel occludes.Because platelet activation by ADP continues in the low shear case, ADP continuesto be released and its concentration remains significant for the entire period leadingup to occlusion.
4. Discussion.
We developed the first mathematical model of flow-mediatedprimary hemostasis in an extravascular injury which is able to track the process frominitial platelet deposition to occlusion. Model calibration with PDE and MFA mod-els yielded excellent agreement in model metrics. ODE model comparisons with thePDE model velocity profiles, amounts of ADP, and total number of bound plateletsvalidate ODE model components describing the fluid, platelets, ADP, and aggregateformation dynamics. The model met the difficult and nuanced challenge of cap-turing injury occlusion in a non-spatial description. We compared MFA occlusiontimes to ODE model output as a function of ADP-independent platelet activationrate ( A ([ ADP ] agg , [ ADP ] BL ) = 1). Next, we studied both ADP-independent andADP-dependent platelet activation and how both activation and flow affect aggregatestructure. Metrics of transport, ADP concentration, and aggregate height and den-sity confirmed that ADP dilution significantly limited aggregate formation. However,we have shown experimentally that platelet aggregates can fully occlude the injurychannel in our bleeding chip (data not shown) in the absence of coagulation, andthus there is a need to extend and further refine our platelet model. In particular,force-dependent bond formation and breaking in the rates of platelet adhesion andcohesion will likely be necessary to include. This will require the modeling of moremechanisms of platelet activation but it will also allow for platelet detachment froma growing aggregate. Most of the model parameters have beeneither previously estimated from the literature or fit to the model output generatedby the PDE and MFA models. The previously estimated parameters are found in
This manuscript is for review purposes only. K.G. LINK, ET AL.
Appendix C (Table 5 and Table 6). The parameter values we are less certain of arenumbers corresponding to the functions describing adhesion and cohesion ( d se , k b ,and s b ) and the permeability constant C K , which prescribes the growing resistance toflow within the aggregate. In reference to d se , k b and s b , the values were determinedthrough tuning ODE model output with flow rates and the number of bound plateletsin the aggregate associated with the PDE model. Details are found in Appendix C.1(Table 7). We determine a value for C K using the protocol described in Appendix C.2.In addition to these constant parameters, the number density of mobile unactivatedplatelets at the inlet to the injury channel, P up ( y ), is prescribed. We assume thatthe shape of the distribution does not evolve as the top and bottom aggregates grow.However, the model framework allows for the change in flow to affect the evolution of P up ( y ), the upstream platelet distribution.In this paper, we consider only one pathway of platelet activation and do notconsider coagulation. In reality, ADP, thrombin (resulting from coagulation), andsubendothelial collagen elicit different types of platelet activation responses. For ex-ample, thrombin cleaves protease-activated receptors on the platelet surface, PAR-1and PAR-4, and is a stronger activator of platelets than ADP [1, 27]. There is in-creasing indication that partially activated and even unactivated platelets may be ableto bind, if only transiently, to a growing aggregate. Such binding would increase aplatelet’s time of exposure to other chemical agonists and possibly lead to more acti-vation and faster aggregate growth. This further motivates the extension of the modelto include alternative mechanisms of platelet activation and binding/unbinding.The processes of platelet binding and unbinding are dictated by the types ofbonds that form between platelets and the vessel wall as well as amongst platelets.Along with the transport of platelets to the site of injury, the flow conditions inthe local environment determine the likelihood that bonds will form and break. Atlow shear rates ( <
200 (1/s)), unactivated platelets near the wall adhere to the col-lagen embedded in the subendothelial matrix [20, 39]. Specifically, bonds betweenGPVI receptors and collagen form and trigger activation of integrin receptors ( α β & α IIb β ). The platelet-collagen bonds bring the platelet to rest and allow for firmadhesion to the subendothelial surface. Adhesion of platelets is a multistep processat high shear rates ( >
500 (1/s)). More specifically, unactivated and mobile plateletstransiently bind to immobilized vWF molecules associated with a thrombus via theirGPIb receptors [19, 35]. This process slows the platelets and allows time for activa-tion of integrins that mediate the formation of strong, long-lived α IIb β -fibrinogenbonds [19, 36, 38]. Future work will include model extensions of this framework thatdifferentiate platelets by the types of bonds that form on their surfaces and by dragon the aggregate that effects individual bonds amongst the aggregated platelets.In addition to extending our occlusive model of platelet aggregation to includemore mechanisms of activation and bond formation and breaking, we can include thebiochemical reactions of coagulation, further maximizing the versatility and reducedcomputational cost of the model. The coagulation cascade occurs on the surfaces ofplatelets producing a key end product thrombin. Thrombin converts fibrinogen tofibrin, which polymerizes forming a stabilizing mesh that supports the platelet aggre-gate. Our previous mathematical model of thrombosis [22, 28, 29] can be linked to theabove model, contributing minimal additional computational cost. Both model exten-sions will utilize the above framework to develop the first flow-mediated mathematicalmodel of hemostasis in an extravascular injury that incorporates both platelet func-tion and coagulation, resulting in a high throughput screening tool that can identifymodifiers of occlusion and embolization. This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. Author Contributions.
Conceptualization:
Kathryn G. Link, Aaron L. Fogelson, Karin Leiderman
Experimental Data:
Matthew G. Sorrells, Keith B. Neeves
Funding acquisition:
Aaron L. Fogelson, Karin Leiderman, Keith B. Neeves
Model Development:
Kathryn G. Link, Aaron L. Fogelson
Validation & Calibration:
Kathryn G. Link, Matthew G. Sorrells, Nicholas A. Danes, KarinLeiderman, Keith B. Neeves, Aaron L. Fogelson
Writing-original draft:
Kathryn G. Link, Aaron L. Fogelson
Writing-review & editing:
Kathryn G. Link, Aaron L. Fogelson, Karin Leiderman, Keith B.Neeves
REFERENCES[1]
H. Andersen, D. L. Greenberg, K. Fujikawa, W. Xu, D. W. Chung, and E. W. Davie , Protease-activated receptor 1 is the primary mediator of thrombin-stimulated platelet pro-coagulant activity , Proceedings of the National Academy of Sciences of the United Statesof America, 96 (1999), pp. 11189–93.[2]
R. Andrews, M. Berndt, and J. Lopez , Platelets , Burlington: Academic, 2001, ch. Theglycoprotein Ib-IX-V complex, pp. 145–163.[3] a. Baurand, P. Raboisson, M. Freund, C. L´eon, J. P. Cazenave, J. J. Bourguignon,and C. Gachet , Inhibition of platelet function by administration of MRS2179, a P2Y1receptor antagonist , European Journal of Pharmacology, 412 (2001), pp. 213–21.[4]
K. Bledzka, S. S. Smyth, and E. F. Plow , Integrin α IIb β : from discovery to efficacioustherapeutic target , Circulation Research, 112 (2013), pp. 1189–1200.[5] T. V. Colace, R. W. Muthard, and S. L. Diamond , Thrombus growth and embolism ontissue factor-bearing collagen surfaces under flow: role of thrombin with and without fibrin ,Arteriosclerosis, thrombosis, and vascular biology, 32 (2012), pp. 1466–1476.[6]
E. L. Cussler , Diffusion Mass Transfer in Fluid Systems , Cambridge University Press, Cam-bridge, Third ed., 2009.[7]
N. Danes and K. Leiderman , A density-dependent FEM-FCT algorithm with application tomodeling platelet aggregation , International Journal for Numerical Methods in BiomedicalEngineering, (2019), p. e3212.[8]
E. C. Eckstein and F. Belgacem , Model of platelet transport in flowing blood with drift anddiffusion terms , Biophysical Journal, 60 (1991), pp. 53–69.[9]
E. C. Eckstein, A. W. Tilles, and F. J. Millero III , Conditions for the occurrence of largenear-wall excesses of small particles during blood flow , Microvascular Research, 36 (1988),pp. 31–39.[10]
P. Elizondo and A. L. Fogelson , A mathematical model of venous thrombosis initiation ,Biophysical Journal, 111 (2016), pp. 2722–2734.[11]
A. Fick , On liquid diffusion , The London, Edinburgh, and Dublin Philosophical Magazine andJournal of Science, 10 (1855), pp. 30–39.[12]
A. L. Fogelson, Y. H. Hussain, and K. Leiderman , Blood Clot Formation under Flow : TheImportance of Factor XI Depends Strongly on Platelet Count , Biophysical Journal, 102(2012), pp. 10–18.[13]
A. L. Fogelson and K. B. Neeves , Fluid mechanics of blood clot formation , Annual Reviewof Fluid Mechanics, 47 (2015), pp. 377–403.[14]
A. L. Fogelson and K. B. Neeves , Fluid Mechanics of Blood Clot Formation , Annual Reviewof Fluid Mechanics, 47 (2015), pp. 377–403.[15]
A. L. Fogelson and N. Tania , Coagulation under flow: The influence of flow-mediated trans-port on the initiation and inhibition of coagulation , Pathophysiology of Haemostasis andThrombosis, 34 (2006), pp. 91–108.[16]
A. R. Gear , Platelet adhesion, shape change, and aggregation: rapid initiation and sig-nal transduction events , Canadian Journal of Physiology and Pharmacology, 72 (1994),pp. 285–294.[17]
E. F. Grabowski, J. T. Franta, and P. Didisheim , Platelet aggregation in flowing bloodin vitro. ii. dependence of aggregate growth rate on adp concentration and shear rate ,Microvascular Research, 16 (1978), pp. 183–195.[18]
B. Hechler and C. Gachet , P2 receptors and platelet function , Purinergic Signalling, 7(2011), pp. 293–303.[19]
S. P. Jackson , Review Article: The growing complexity of platelet aggregation , Blood, 109(2007), pp. 5087–5096.[20]
S. P. Jackson, W. S. Nesbitt, and S. Kulkarni , Signaling events underlying thrombusformation , Journal of Thrombosis and Haemostasis, 1 (2003), pp. 1602–1612.
This manuscript is for review purposes only. K.G. LINK, ET AL.[21]
C. KJ and C. Clemetson , Platelets , Burlington: Academic, 2007, ch. Platelet Receptors,pp. 201–213.[22]
A. L. Kuharsky and A. L. Fogelson , Surface-Mediated Control of Blood Coagulation : TheRole of Binding Site Densities and Platelet Deposition , Biophysical Journal, 80 (2001),pp. 1050–1074.[23]
K. Leiderman, W. C. Chang, M. Ovanesov, and A. L. Fogelson , Synergy Between TissueFactor and Exogenous Factor XIa in Initiating Coagulation , Arteriosclerosis, Thrombosis,and Vascular Biology, 36 (2016), pp. 2334–2345.[24]
K. Leiderman and A. L. Fogelson , Grow with the flow: A spatial-temporal model of plateletdeposition and blood coagulation under flow , Mathematical Medicine and Biology, 28(2011), pp. 47–84.[25]
K. Leiderman and A. L. Fogelson , The influence of hindered transport on the developmentof platelet thrombi under flow , Bulletin of mathematical biology, 75 (2013), pp. 1255–1283.[26]
K. Leiderman and A. L. Fogelson , The Influence of Hindered Transport on the Developmentof Platelet Thrombi Under Flow , Bulletin of Mathematical Biology, 75 (2013), pp. 1255–1283.[27]
Z. Li, M. K. Delaney, K. A. O’Brien, and X. Du , Signaling during platelet adhesion andactivation , Arteriosclerosis, Thrombosis, and Vascular Biology, 30 (2010), pp. 2341–2349.[28]
K. G. Link, M. T. Stobb, J. Di Paola, K. B. Neeves, A. L. Fogelson, S. S. Sindi, andK. Leiderman , A local and global sensitivity analysis of a mathematical model of coagu-lation and platelet deposition under flow , PLOS ONE, 13 (2018).[29]
K. G. Link, M. T. Stobb, M. G. Sorrells, M. Bortot, K. Ruegg, M. J. Manco-Johnson,J. A. Di Paola, S. S. Sindi, A. L. Fogelson, K. Leiderman, et al. , A mathematicalmodel of coagulation under flow identifies factor v as a modifier of thrombin generationin hemophilia a , Journal of Thrombosis and Haemostasis, (2019).[30]
W. McCabe, J. Smith, and P. Harriott , Unit Operations of Chemical Engineering , CivilEngineering, McGraw-Hill Education, 2005.[31]
R. W. Muthard and S. L. Diamond , Blood clots are rapidly assembled hemodynamic sensors:flow arrest triggers intraluminal thrombus contraction , Arteriosclerosis, Thrombosis, andVascular Biology, 32 (2012), pp. 2938–2945.[32]
K. Nogami and M. Shima , Phenotypic heterogeneity of hemostasis in severe hemophilia , inSeminars in thrombosis and hemostasis, vol. 41, Thieme Medical Publishers, 2015, pp. 826–831.[33]
P. Ohlmann, S. de Castro, G. G. Brown, C. Gachet, K. A. Jacobson, and T. K.Harden , Quantification of recombinant and platelet P2Y1receptors utilizing a [125I]-labeled high-affinity antagonist 2-iodo-N6-methyl-(N)-methanocarba-2’-deoxyadenosine-3,5-bisphosphate ([125I]MRS2500) , Pharmacological Research, 62 (2010), pp. 344–351.[34]
K. Rana and K. B. Neeves , Blood flow and mass transfer regulation of coagulation , BloodReviews, 30 (2016), pp. 357–368.[35]
Z. M. Ruggeri , Platelet adhesion under flow. , Microcirculation (New York, N.Y. : 1994), 16(2009), pp. 58–83.[36]
Z. M. Ruggeri, J. a. Dent, and E. Sald´ıvar , Contribution of distinct adhesive interactionsto platelet aggregation in flowing blood. , Blood, 94 (1999), pp. 172–8.[37]
K. S. S. Kobayashi and D. N. Ku , Occlusive thrombus growth at high shear rates: comparisonof whole blood and platelet rich plasma at constant pressure , J. Thromb. Haemost., 14(2016), pp. 6–7.[38]
B. Savage, F. Almus-Jacobs, and Z. M. Ruggeri , Specific synergy of multiple substrate-receptor interactions in platelet thrombus formation under flow , Cell, 94 (1998), pp. 657–666.[39]
B. Savage, E. Sald´ıvar, and Z. M. Ruggeri , Initiation of platelet adhesion by arrest ontofibrinogen or translocation on von willebrand factor , Cell, 84 (1996), pp. 289–297.[40]
R. Schoeman, K. Rana, N. Danes, M. Lehmann, J. Di Paola, A. Fogelson, K. Leiderman,and K. Neeves , A microfluidic model of hemostasis sensitive to platelet function andcoagulation , Cellular and Molecular Bioengineering, 10 (2017), pp. 3–15.[41]
L. F. Shampine and M. W. Reichelt , The matlab ode suite , SIAM Journal on ScientificComputing, 18 (1997), pp. 1–22.[42]
T. Skorczewski, L. C. Erickson, and A. L. Fogelson , Platelet motion near a vessel wallor thrombus surface in two-dimensional whole blood simulations , Biophysical journal, 104(2013), pp. 1764–1772.[43]
G. Tangelder, H. C. Teirlinck, D. W. Slaaf, and R. S. Reneman , Distribution of bloodplatelets flowing in arterioles , American Journal of Physiology-Heart and Circulatory Phys-iology, 248 (1985), pp. H318–H323.
This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. [44] A. W. Tilles and E. C. Eckstein , The near-wall excess of platelet-sized particles in bloodflow: its dependence on hematocrit and wall shear rate , Microvascular Research, 33 (1987),pp. 211–223.[45]
L. Timmermans, P. Minev, and F. Van De Vosse , An approximate projection scheme forincompressible flow using spectral elements , International journal for Numerical Methodsin Fluids, 22 (1996), pp. 673–688.[46]
V. T. Turitto and E. F. Leonard , Platelet adhesion to a spinning surface , ASAIO Journal,18 (1972), pp. 348–354.[47]
L. Valentino , Blood-induced joint disease: the pathophysiology of hemophilic arthropathy ,Journal of Thrombosis and Haemostasis, 8 (2010), pp. 1895–1902.[48]
H. H. Versteeg, J. W. M. Heemskerk, M. Levi, and P. H. Reitsma , New Fundamentals inCoagulation , Physiology Review, 93 (2013), pp. 327–358.[49]
A. Wufsus, N. Macera, and K. Neeves , The hydraulic permeability of blood clots as afunction of fibrin and platelet density , Biophysical Journal, 104 (2013), pp. 1812–1823.[50]
C. Yeh, A. C. Calvez, and E. C. Eckstein , An estimated shape function for drift in aplatelet-transport model , Biophysical Journal, 67 (1994), pp. 1252–1259.
Appendix A. Flow Through Bleeding Chip.A.1. Circuit Analog System.
We use hydrodynamic analogs of Kirchoff’s andOhm’s laws to formulate a linear system of equations that describes the flow throughthe hydraulic circuit (HC) shown in Figure 3A. Q − Q M = Q , (A.1) Q − Q M = Q , (A.2) P B − P = Q R , (A.3) P W − P = Q R , (A.4) P − P CB = Q R , (A.5) P − P CW = Q R , (A.6) P − P = Q M R M . (A.7)The inlet and outlet pressures P B , P W , P CB , and P CW are known and we assume theresistance through the injury channel R M is an input parameter determined by theBrinkman-Stokes calculation. Additional input parameters include R , R , R , and R , which correspond to the resistances associated with the blood and wash channels.Equations (A.1)-(A.7) yields a system of seven equations with seven unknowns. Wecan write them as the following matrix system.(A.8) R R R R R M Q Q Q Q Q M P P = P B P W − P CB − P CW Values of pressures and resistances used in the hydraulic circuit (HC) system are sum-marized in Table 1, Table 2, and Table 3. There are variable viscosities correspondingto the blood, wash and injury channels. Furthermore, the downstream end of the washchannel will be a mix of the blood and the wash when blood can flow freely from theblood channel through the injury channel. We account for these viscosities throughthe resistance terms found in Table 2. Assume that the viscosity of the downstreamend of the wash channel is an average of the blood and wash viscosities. In the fluidsimulations, we use a volume-averaged viscosity through the entire domain µ = 0 . This manuscript is for review purposes only. K.G. LINK, ET AL. dynes · s/cm . We also assume there is an effective length of the injury channel thatspans outside its 150 µ m length, corresponding to the laminar flow that enters andexits the injury channel. The effective length of the channel, L eff h = 165 µ m, is usedin the calculation of the pressure gradient through the injury channel when using thePDE ‘flow map’ output. Table 1
Dimensions of Hydraulic Circuit (HC).
See Figure A-B for reference.
Length from P B to P (and P W to P ) 85 µ mLength from P to P CB (and P to P CW ) 85 µ mWidth of vertical channels 100 µ mLength of injury channel 150 µ mWidth of injury channel 20 µ mDepth of injury channel 60 µ m Table 2
HC Viscosities and Computed Resistances.
See Figure A for reference.
Blood Channel Viscosity 3 . × − PoiseUpstream Wash Channel Viscosity 1 × − PoiseInjury Channel Viscosity 2 . × − PoiseDownstream Wash Channel Viscosity 2 . × − Poise R , R . × Pa s/m R M (no aggregate) 6 . × Pa s/m R . × Pa s/m R . × Pa s/m Table 3
Pressures used in the HC.
See Figure A for reference.
Pressure Inlet (Blood, P B ) 6 . × PaPressure Inlet (Wash, P W ) 1 . × PaPressure Outlet (Blood, P CB ) 2 . × PaPressure Outlet (Wash, P CW ) 0 PaPressure Drop Across Injury P − P (no aggregate) 2 . × Pa Table 4
Parameters that characterize the geometry and flow, mass transfer, and reactionregimes in the injury channel.
The average blood velocity U , ν kinematic viscosity, D diffusivityof platelets, γ wall shear rate. ** D h , hydraulic diameter [2 d h L z / ( d h + L z )] . Parameters Expression Value Constants
Aspect ratio d h /L z . d h = 20 µ m, L z = 60 µ mRelative injury size L h /d h . L h = 150 µ mReynolds number ( R e ) Ud h /ν . ν = 1 . × − m / s, U = 2 . × − m/sEntrance Length ( L e ) 0 . R e D h . µ m D h = 30 µ mPeclet number ( P e ) γd h / D p . × γ = 2400 s − , D p =2 . × − m /sMass transfer coeff. ( k m ) ( D ρ/ L h ) / . × − m/s ρ = 1060 kg/m A.2. Experimental Methods & Materials.A.2.1. Materials.
Bovine serum albumin (BSA) (A9418), 3,3’-dihexylox acar-bocyanine iodide (DiOC6) (D273) were from Sigma-Aldrich (St Louis, MO, USA).3.2% sodium citrate vacutainers (369714) and 21-gauge Vacutainer ®
21 Safety-Lokblood collection sets (367281) were from Becton Dickson (Frankwood Lakes, NJ,
This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY.
A.2.2. Blood Collection and Preparation.
Blood was collected from healthydonors by venipuncture into 4.5 mL PPACK vacutainers (final concentration of 75 µ M) using a 21 gauge needle. Prior to collecting blood into the PPACK vacutainers,one vacutainer of blood was collected into a 3.2% sodium citrate vacutainer anddiscarded. DiOC6 was added to blood to a concentration of 1 µ M and was incubated at37 ◦ C for 10 min prior to the microfluidic device assay. The study and consent processwas approved from the Colorado Multiple Institutional Review Board in accordancewith the Declaration of Helsinki.
A.2.3. Device Fabrication, Preparation, and Operation.
A PDMS mi-crofluidic device was created using stand photolithography and soft lithography tech-niques. The device was made with an ‘H-shaped’ geometry, with a ‘blood’ chan-nel (w=110 µ m, h=50 µ m) on the right, a ‘wash’ channel on the left (w=110 µ m,h=50 µ m), and an ‘injury’ channel (l=150 µ m, w=20 µ m, h=50 µ m) connecting thetwo vertical channels (Figure 1A). The device was plasma bonded to a PDMS coatedglass slide, and a solution of the collagen related peptides CRP, GFOGER, and VWF-BP were patterned in the injury channel and wash channel at a concentration of 250 µ g/mL for each peptide for 4 hours. Following patterning, the entire device was rinsedand blocked with 2% BSA in PBS for 45 min. A.2.4. Whole Blood Microfluidic Flow Assays.
Pressure controllers (Flu-igent MFCS-EZ) coupled to flow meters (Fluigent FRP) were used to perfuse bloodand wash buffer at a flowrate of 5.5 µ L/min and 10 µ L/min, respectively, through thedevice. Images of thrombus in the injury channel were obtained using brightfield andepifluorescent microscopy on an inverted microscope (Olympus IX81, 20X Objective,NA 0.4, Hammamatsu Orca R2 Camera). Assays were run until the injury channelwas fully occluded. The time to occlusion in the device was determined through avariety of ways. Optically, we defined the time to occlusion as the time point in whichno red cells could be seen passing through the injury channel. Additionally, we usedthe flow meter measurements on the inlet and outlet of the device to compute theflow rate through the injury channel. This flow measurement was used to calculatetwo other occlusion metrics: the time to reach 10% of the endpoint flow rate and thetime needed to reach 0.35 µ L/min flow.
A.3. Creation of flow map.
For comparisons with the PDE model, the ODEmodel utilizes a flow rate through the injury channel to prescribe the appropriate cor-responding pressure drop and this is needed for all possible combinations of aggregateheights and densities. To acquire these approximations, we used a two-dimensionalPDE model of flow through an extravascular injury channel that was initialized tobe filled with various heights and densities of porous material to represent possibleaggregate formations. For each of the various aggregates, the model was simulated
This manuscript is for review purposes only. K.G. LINK, ET AL. until the flow came to an equilibrium state and resulting flow rates were recorded; theresults were assimilated into what we call the ‘flow map’. Details of these calculationsare as follows.The geometric configuration for the 2D model is the ‘H’ domain as depicted inFigure 1A and as described in detail in our previous work [7]. Briefly, there aretwo vertical channels through which blood and wash (buffer) flow. These channelsare bifurcated by a horizontal injury channel and, with specified flow rates throughthe vertical channels, blood flows through the injury channel and out into the washchannel. Within the injury channel, porous aggregates can be placed and fixed tostudy pressure drops and flows, or the aggregates can grow dynamically accordingto platelet and ADP equations. Fluid at the fixed flow rate enters the domain atthe top of the vertical channels and conditions on pressure at the bottom of thevertical channels are prescribed. For studies with fixed, initial distributions of porousaggregate, the aggregate is considered only on the interior of the horizontal injurychannel and is placed in such a way that the nonzero-porosity is centered along alength L h − x , where ∆ x is the size of one computational grid cell. For thesesimulations, ∆ x ≈ µ m.Fluid throughout the computational domain is modeled with the Navier-Stokes-Brinkman equations to account for unimpeded flow in the vertical channels and im-peded flow through the porous aggregate within the injury channel. The equationstake the form: ρ (cid:18) ∂(cid:126)u∂t + (cid:126)u · ∇ (cid:126)u (cid:19) = −∇ p − µ ∇ (cid:126)u − µα ( θ B ) (cid:126)u, (A.9) ∇ · (cid:126)u = 0 , (A.10)where (cid:126)u and p represent fluid velocity and pressure, respectively, µ is the dynamicviscosity and ρ is the fluid density. The term − µα ( θ B ) (cid:126)u represents a frictional resis-tance to the fluid motion due to the presence of the platelet aggregate with volumefraction of bound platelets, θ B . Equations (A.9)-(A.10) are solved numerically usingthe rotational projection method, as described previously [7, 45]. We assume a func-tional form of the friction term, α ( θ B ), that follows the Carman-Kozeny relation andthe particular value of C K is specified in Appendix C.2.Flow rates that correspond to experimental setups at the top of the blood channel Q b and wash channel Q w are used to set the velocity boundary conditions at theirrespective inlets by using a parabolic velocity profile of the form: u y = 32 Q i W ( x − w i ) W (cid:18) − x − w i W (cid:19) where u y is the y − component (parallel to the vertical channels) of the velocity field (cid:126)u , W is the width of the blood channel, w i the center point of inlet channel correspondingto index i , and Q i is the given flow rate corresponding to index i , with i = b, w . Outletpressure conditions are set using an open boundary condition of the form: ∂u n ∂n − p = s, ∂u t ∂n = 0where u n and u t are the normal and tangential components of the fluid velocity,respectively, p is the fluid pressure, and s is the prescribed (constant) pressure givenfrom the circuit calculation, as described above. This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. d b , d t ∈ { , , , , µm } ,where d b , d t are the height of the bottom and top aggregate, respectively, and for θ B b , θ B t ∈ [0 . , . , · · · . θ B b , θ B t are the density of the bottom and top ag-gregate, respectively. The resulting flow map acts as a “look-up” table; given ODEmodel values of d b , d t , θ B b , θ B t we do quadrilinear interpolation from table values todetermine Q P DE ( d b , d t , θ B b , θ B t ). A.4. PDE Model of Platelet Aggregation.
In this section, we detail all thecomponents of the PDE model of platelet aggregation developed in tandem with theODE model described in this paper. The evolution equations for the three species ofplatelets are given by: ∂P m,u ∂t = −∇ · (cid:2) W ( θ T ) ( (cid:126)uP m,u − D p ∇ P m,u ) (cid:3)(cid:124) (cid:123)(cid:122) (cid:125) advection and diffusion − k adh( x ) (cid:0) P max − P b,a (cid:1) P m,u (cid:124) (cid:123)(cid:122) (cid:125) adhesion − k act ([ADP]) P m,u (cid:124) (cid:123)(cid:122) (cid:125) activation by ADP , (A.11) ∂P m,a ∂t = −∇ · (cid:2) W ( θ T ) ( (cid:126)uP m,a − D p ∇ P m,a ) (cid:3) − k adh ( x ) (cid:0) P max − P b,a (cid:1) P m,a − k coh η ( x ) P max P m , a (cid:124) (cid:123)(cid:122) (cid:125) cohesion + k act ([ADP]) P m,u , (A.12) ∂P b,a ∂t = k adh ( x ) (cid:0) P max − P b,a (cid:1) P m,a + k coh η ( x ) P max P m , a , (A.13)where θ T = (cid:0) P m,u + P m,a + P b,a P max (cid:1) θ max , and D p , k adh ( x ), k coh , and W ( θ T ) are the plateletdiffusion coefficient, spatially-dependent adhesion rate, cohesion rate coefficient, andhindered platelet flux coefficient, all used as previously described [7, 24]. The pa-rameter, η ( x ) is a local approximation to the bound platelet fraction that indicatesboth the location and density of bound platelets to mobile platelets as part of thecohesion rate; the cohesion rate is enhanced where the local bound platelet fractionis increased. We approximate this parameter on a triangular but structured mesh,taking advantage of the interpolated values of θ B , which we will denote as ˜ θ B , definedusing the finite element method. Since θ B is defined on the finite element mesh withpiecewise linear interpolants, we can approximate θ B anywhere in the domain by sim-ply evaluating these interpolants. We allow η ( x ) at nodes to depend on values of θ B that are up to h p = 2 µ m microns away. Given a nodal point ( x j , y j ) on the triangularmesh, we use a weighted average that depends the most strongly on θ B ( x j , y j ): This manuscript is for review purposes only. K.G. LINK, ET AL. η ( x j , y j ) = 14 θ B ( x j , y j ) + 18 (cid:104) ˜ θ B ( x j + h p , y j ) + ˜ θ B ( x j − h p , y j )+ ˜ θ B ( x j , y j + h p ) + ˜ θ B ( x j , y j − h p ) (cid:105) , + 116 (cid:104) ˜ θ B ( x j + h p , y j + h p ) + ˜ θ B ( x j − h p , y j + h p )+ ˜ θ B ( x j + h p , y j − h p ) + ˜ θ B ( x j − h p , y j − h p ) (cid:105) . Tracking ADP concentration and its effects on platelet activation in the PDE modelis similar to that in the ODE model. Activation occurs through(A.14) k act = k adp [ ADP ][ ADP ] + [
ADP ] ∗ where k adp = 3 . − is the activation rate and [ ADP ] ∗ is the critical concentration ofADP for which there is significant activation of platelets by ADP. ADP molecules arereleased by bound, activated platelets and are transported in the fluid by advectionand diffusion. Hence, we track the ADP concentration with an advection-diffusion-reaction equation:(A.15) ∂ [ADP] ∂t = −∇ · ( (cid:126)u [ADP] − D a ∇ [ADP]) + σ release where D a is the diffusion coefficient for ADP, and σ release is the release rate, and theyare defined exactly the same way as is done in the ODE model, described above. Thecoupled flow equations are the Navier-Stokes-Brinkman equations described abovewith the volume fraction due to bound platelets defined as θ B = (cid:0) P b,a P max (cid:1) θ max . Appendix B. Stokes-Brinkman Velocity System.
To determine the veloc-ity profile through the injury channel described in Subsection 2.1.2, we solve a linearsystem of six equations and six unknowns. Let y t = d h − d t . A b + B b = − G h (cid:18) µα b (cid:19) , (B.1) A t exp (cid:0) √ α t d h (cid:1) + B t exp (cid:0) − √ α t d h (cid:1) = − G h (cid:18) µα t (cid:19) , (B.2) Ad b + B − A b exp (cid:0) √ α b d b (cid:1) − B b exp (cid:0) − √ α b d (cid:1) = G h (cid:18) µα b + d µ (cid:19) , (B.3) Ay t + B − A t exp (cid:0) √ α t y t (cid:1) − B t exp (cid:0) − √ α t y t (cid:1) = G h (cid:18) µα t + y µ (cid:19) , (B.4) A − A b √ α b exp (cid:0) √ α b d b (cid:1) + B b √ α b exp (cid:0) − √ α b d b (cid:1) = G h (cid:18) d b µ (cid:19) , (B.5) A − A t √ α t exp (cid:0) √ α t y t (cid:1) + B t √ α t exp (cid:0) √ α t y t (cid:1) = G h (cid:18) y t µ (cid:19) . (B.6)The pressure gradient G h , and α b , α t , µ , d b , d t , and the width of the injury channel, d h are prescribed. The unknowns are A , B , A b , B b , A t , B t . This system (B.1)-(B.6) This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. C ABA b B b A t B t = G h µ - α b - α t α b + d α t + y d b y t where the matrix C is defined as C = e √ α t d h e −√ α t d h d b e √ α b d b - e −√ α b d b y t e √ α t y t - e −√ α t y t √ α b e √ α b d b √ α b e −√ α b d b √ α t e √ α t y t √ α t e −√ α t y t . Note that the right-hand side of the system is linearly dependent on the pressuregradient G h . Therefore, we assume G h = 1 and solve the linear system of six equationsfor the six constants that prescribe the flow velocity profile. The velocity profile isused to calculate the resistance R M through the injury channel as described below R M = L h L z (cid:82) d h u ( y ) dy . The effects of the growing aggregate on the flow velocity through the injury channelare captured in the quantity R M , which is an input parameter in the hydraulic circuit(HC) system describing the flow of blood through the bleeding chip. The true G h isdetermined by solving the HC system and is used to scale the velocity profile thatwas obtained with G h = 1. Appendix C. Model Parameters.
Table 5
Diffusion coefficients, ADP granule release, and upstream platelet number density.
Species Value ReferencePlatelets 2 . × − cm /s [46]ADP 5 × − cm /s [17]ˆ A − moles of ADP/plt [24] P up, ∗ . plt/ µ L [24]
C.1. Calibrated Parameters.
Both adhesion and cohesion parameters de-tailed in Subsection 2.2.2 were tuned with the PDE and MFA models. Additionally,constants c and c used in the definitions of the diffusive fluxes of ADP described inSubsection 2.3.3 were determined through calibration with the PDE model. Valuesare summarized in Table 7. C.2. Aggregate Permeability.
When considering aggregate permeability, thesignificant differences between ‘white’ and ‘red’ clots must be considered. A white clotis usually formed under high shear and is predominantly made up of platelets whereas
This manuscript is for review purposes only. K.G. LINK, ET AL.
Table 6
Platelet transitions.
Transition Initial State Final State M − s − s − Unactivated plateletadhering to SE P m,u P b,a k adh = 2 . × [22, 24]Activated plateletadhering to SE P m,a P b,a k adh = 2 . × Activated plateletcohering to boundplatelet P m,a P b,a k coh × P max =1 × [24]Platelet activationby ADP P m,u P m,a k adp = 3 . Table 7
Calibrated parameters.Description Parameter Value
Partition constant k b s b d se µ mGradient constants c , c /
2, 1 / a fibrin-rich red clot forms under low shear or static conditions. Table 8 lists the com-puted permeability comparison results from the literature. It is shown that white clotsare significantly more permeable than red clots. In addition to the type of clot, theflow conditions and experimental setup generate significant variability. For exampleresults in [49] were generated by measuring the permeability of platelet-rich aggregatesat various volume fractions. Specifically, the authors induced clotting with thrombinin suspension with known platelet densities and then measured the permeation of abuffer solution through the clots at defined pressure gradients (constant pressure).Permeability values of 1 . × − ± . × − − . × − ± . × − µ m were cal-culated from platelet-rich clots (PRC) with platelet volume fractions of 0 . − .
61. Inconclusion, both flow conditions and experimental set up yield a wide range of perme-abilities. Given the variability in permeability data found in the literature, we tuned
Table 8
Permeability comparison.Clot Type Permeability [ µ m ] Reference Red 5 . × − Muthard and Diamond (2012) [5]Red 1 . × − Wufsus et al. (2013) Platelet 61% [49]Red 6 . × − Wufsus et al. (2013) Platelet 31% [49]White 0 . . the constant C K using input from both PDE and MFA models. More specifically, weassume a functional form of the Brinkman coefficient that follows the Carman-Kozenyrelation [30] α ( θ B ) = C K ( θ B ) (1 − θ B ) , C K = 3 . × cm − . This particular value of C K was determined by performing flow simulations withthe PDE model in which the injury channel was completely filled with aggregateat a minimum porosity until it matched the experimentally observed flow rate of This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. . µ L/min.
Appendix D. Non-uniform Aggregate Formation.
We considered uniform,near-wall, and shifted upstream platelet distributions P up ( y ) found in Figure 12Aand compared the formation of the resulting aggregates. Additionally, we examinedthe effects of an asymmetric upstream platelet distribution shown in Figure 12B. Thedistributions are defined by the distance y from the bottom wall of the injury channel.Therefore, y = 20 µ m is the location associated with the top wall. Both sets of in silico experiments utilized ADP-independent activation where A ([ ADP ] agg , [ ADP ] BL ) = 1. (A) (B) y ( µ m) P l t . C o n c . ( × / µ L ) UniformNear-wallShifted y ( µ m) P l t . C o n c . ( × / µ L ) Fig. 12 . Upstream Platelet Distributions: (A)Uniform ( P up ( y ) = P up, ∗ ), near-wall, andshifted upstream platelet distributions. These are motivated by in vivo [43] and in vitro [8, 9, 44, 50]observations. Additional details are found in [24]. Under both high and low shear conditions, aggregates simulated with the uniformplatelet distribution are the slowest growing in both thickness and density (Figure 13A-B). The use of a near-wall platelet distribution, which has the highest concentrationof platelets closest to the walls, results in the fastest occlusion time of approximately1.75 minutes. The shifted platelet distribution produced aggregates that occlude ap-proximately 30-45 seconds slower. A similar trend is seen in the low shear conditionsFigure 13 C-D, however occlusion times are longer. The results from these symmetriccases yield aggregates with varying thicknesses and densities, identifying the upstreamplatelet distribution P up ( y ) as a modifier of aggregate formation. Additionally, wewere motivated to perform experiments that address the effects of asymmetric up-stream platelet distribution on aggregate formation.Figure 14 A-D show aggregate thicknesses and densities associated with the bot-tom (purple) and top (green) walls. We chose an asymmetric distribution found inFigure 12 B, where there is a greater concentration of platelets near the bottom wallof the injury channel. The aggregate grows asymmetrically in both thickness anddensity. In high shear (HS) conditions, the bottom aggregate region is thicker anddenser after 1.25 minutes. There is, however, a significant increase in the densityof the top aggregate after 1.25 minutes, resulting in a thin but dense top aggregate.Similar behavior is seen in results under low shear (LS) conditions, however, thereis a delay. Therefore, the model is able to produce aggregates of different orienta-tions, thicknesses and densities by varying model inputs such as upstream plateletdistribution and flow conditions. This manuscript is for review purposes only. K.G. LINK, ET AL. (A) (B)
Time (min.) d agg ( m i c r on s ) Aggregate Thickness - HS
UniformNear WallShifted
Time (min.) B Agggregate Density - HS
UniformNear WallShifted (C) (D)
Time (min.) d agg ( m i c r on s ) Aggregate Thickness - LS
UniformNear WallShifted
Time (min.) d agg ( m i c r on s ) Aggregate Thickness - LS
UniformNear WallShifted
Time (min.) B Aggregate Density - LS
UniformNear WallShifted
Time (min.) B Agggregate Density - LS
UniformNear WallShifted
Fig. 13 . Symmetric Aggregate Formation:
Given different symmetric upstream plateletdistributions P up ( y ) : uniform, near-wall, and shifted (see Figure A) with ADP-independentplatelet activation A ([ ADP ] agg , [ ADP ] BL ) = 1 under high shear (HS) conditions ( γ = 2400 s − ) A) the aggregate thickness B) and density vary as functions of time. Similar comparisons underlow shear (LS) conditions ( γ = 100 s − ) yield aggregates of varying thickness (C) and density (D) .Insets show results after 20 minutes of simulation time confirming occlusion. This manuscript is for review purposes only.
LATELET AGGREGATION IN AN EXTRAVASCULAR INJURY. (A) (B) Time (min.) d agg ( m i c r on s ) Aggregate Thickness - HS
BottomTop
Time (min.) B Aggregate Density - HS
BottomTop (C) (D)
Time (min.) d agg ( m i c r on s ) Aggregate Thickness - LS
BottomTop
Time (min.) B Aggregate Density - LS
BottomTop
Fig. 14 . Asymmetric Aggregate Formation:
Given an asymmetric upstream platelet dis-tribution (see Figure B) with ADP-independent platelet activation ( A ([ ADP ] agg , [ ADP ] BL ) = 1 )and high shear conditions ( γ = 2400 − ), the aggregate thickness A) and aggregate density B) asfunctions of time. Similar comparisons under low shear (LS) conditions ( γ = 100 s − ) yield aggre-gates of varying thickness (C) and density (D) ..