Tissue failure propagation as mediated by circulatory flow
Gurdip Uppal, Gokhan Bahcecioglu, Pinar Zorlutuna, Dervis Can Vural
TTissue failure propagation as mediated by circulatory flow
Gurdip Uppal , ∗ , Gokhan Bahcecioglu , ∗ , Pinar Zorlutuna , ∗∗ , Dervis Can Vural , ∗∗ Department of Physics, University of Notre Dame, IN, USA Department of Aerospace and Mechanical Engineering, University of Notre Dame, IN, USA ∗ Equally contributing first authors ∗∗ Equally contributing corresponding authors ∗∗ Correspondence: [email protected], [email protected]
Aging is driven by subcellular processes that are relatively well-understood. However the qualitative mech-anisms and quantitative dynamics of how these micro-level failures cascade to a macro-level catastrophe ina tissue or organs remain largely unexplored. Here we experimentally and theoretically study how cell failurepropagates in a synthetic tissue in the presence of advective flow. We argue that cells secrete cooperative factors,thereby forming a network of interdependence governed by diffusion and flow, which fails with a propagatingfront parallel to advective circulation.
Significance: Mortality rates typically increase for complex organisms as they age. This leads us to suggest that agingdepends on interactions between cells. As cells become damaged, the effect propagates to other cells, eventually leadingto a systemic catastrophe. Yet it is unclear how this failure dynamically propagates. Here we present experimentswith synthetic tissues and analogous analytical models to investigate the dynamics of failure propagation. Our maincontribution is a detailed investigation of failure propagation when interactions are mediated by advective flow. Wefind analytical expressions for when a pronounced propagation occurs, its velocity, and acceleration in terms of systemparameters.
INTRODUCTION
The aging and death of an organism is typically attributed tosubcellular mechanisms such as reactive oxygen species dam-age or slowing down of tissue repair due to shortening telom-eres. However, an organism does not die because it graduallyruns out of cells, but rather, cellular level failures cascade totissues and organs that lead to a relatively sudden systemiccatastrophe. While microscopic mechanisms of cellular mal-function are relatively well studied [1–4], how failure spreadsfrom the subcellular level to tissues and organs and ultimatelythe organism is largely unknown.In [5] the failure of an organism was modeled as a reliabil-ity circuit where cells within an organ are connected by ORgates, so that an organ fails when all its cells fail), and theorgans are connected by AND gates, so that organism diesas soon as one organ fails. In [6], aging was described asfailures taking place on a complex network of interdependentbuilding blocks. Here, when a node in the network malfunc-tions, so will those that depend on it. As a result, few ran-dom microscopic failures can propagate into many others, ul-timately leading to a catastrophe. [6] could bridge micro scalemalfunctions with their experimentally observed macroscopicmanifestations such as organismic death and population de-mographics, and fit experimental data such as [7] and [8] (alsocf. appendix of [9]). The network models have also been in-sightful in studying frailty [10, 11].While describing a complex organ or an entire organismas a random network of interdependencies is a useful startingpoint to understand how it fails [12], it is also a rather crudeoversimplification. First, in real biological systems, the large-scale structure of interdependence network is far from “ran-dom”. Secondly, there can be varying amounts and varying types of dependencies. In an actual complex biological sys-tem, the exchange of signals and goods between cells occursvia either diffusion or complex patterns of vascular circula-tion. As such, biophysically grounded analogs of interdepen-dence networks are necessary.Such biophysical extensions have been investigated experi-mentally [13] and theoretically [14] to understand how failurepropagates through tissues, as mediated by the loss of diffus-ing cooperative factors. These cooperative factors could becytokines (e.g. interleukin 15), growth factors (e.g. epidermalgrowth factor), survival factors (e.g. insulin-like growth fac-tor 1), and antioxidant enzymes (e.g. superoxide dismutase3) [15–23] diffusing across cells. However, the role of con-vective circulation, which is the primary mode of transport inlarge complex organisms, is missing.The central argument of this paper is that cell failure cas-cades to higher structures along the circulatory network of anorganism, in the same direction as the advective flow. Tothis end, we theoretically and experimentally study the me-chanical nature of damage propagation through tissues in thepresence of advective flow, and demonstrate and quantify howfailures accumulate and propagate in relation to the flow direc-tion. While our experiments are conducted on synthetic mam-malian tissues in a microfluidic device, from here and withthe help analytical arguments, we aim to derive more generallessons about the aging dynamics in complex organisms invivo.The demographic hallmark of aging is a monotonically in-creasing mortality rate µ ( t ) . In strongly aging species (likehumans), the probability of death of an old individual is manytimes larger than a young individual. Interestingly, thereare phylogenetic correlations in aging characteristics [24]:Mammalian populations have the steepest mortality curves, a r X i v : . [ q - b i o . T O ] M a y whereas µ ( t ) of amphibians and reptiles change little overtime, and plants tend to age even less, and they can exhibitmortality rates that even decrease over time. The cause be-hind these phylogenetic trends is not entirely clear, but onepossibility may be the differences in how goods and signalsare transported between cells, manifesting as differences inhow failures accumulate and propagate.For example, assaulting a sizeable portion of a shrub mightnot kill it, because of how little the rest of the shrub dependson the assaulted portion. The interdependence structure of ashrub is highly localized. In organisms with fast and efficientcirculatory systems, however, advective transport enables anycell to depend on any other, no matter how distant. A mal-function in an animal gland or organ will affect all cells thatare coupled to it via bloodstream. Thus, as much as convec-tive flow propagates goods and signals, it should also propa-gate failure along the same path once the goods and signals gomissing.In our experiments, we encapsulate rat cells in a non-degradable hydrogel (where they cannot proliferate, migrate,or contact each other), seated in an engineered microfluidicmicrochannel, through which we apply unidirectional flow ofmedia, meant to emulate blood, lymph or interstitial flow. Wethen analyze cell death rate across the channel as a function oftime and flow rate. Directing cooperative factors downstreamleads to a wave of failure starting near the inlet and propa-gating towards the outlet, which we analytically describe us-ing a mechanistic diffusion-convection-population-dynamicsmodel.From this model we obtain relevant length scales in termsof system parameters and determine the conditions for wherea failure wave should originate, and obtain the velocity andacceleration of its propagation. MATERIALS AND METHODSConstruction of the tissue engineered model
Fabrication of the microfluidic device.
Poly(dimethylsiloxane) (PDMS) prepolymer and curingagent (Dow Corning, USA) were mixed thoroughly at a 10:1ratio and cast on a silicon-base mold with 20mm × × × width × height) ridges, and cured at 80°Cfor 30min. The PDMS and a glass slide were treated with airplasma for 1min and immediately bound to each other withthe channel side facing the slide. The device was sterilizedunder UV for h. Preparation of the polyethylene glycol (PEG) solu-tion.
Maleimide polyethylene glycol (PEG) succinimidyl car-boxymethyl ester (PEG-NHS, Mw 3400Da, JenKem Technol-ogy) was conjugated with tyrosine-arginine-glycine-asparticacid-serine (Y
RGD
S, Bachem) as described previously, to ob-tain the PEG-RGD [13]. The PEG-RGD and 4-arm PEG-acrylate (4-PEG-ACR, 20kDa, JenKem Technology) weremixed at a 1.5:8.5 (w/w) ratio, and a 20% (w/v) solutionwas prepared in culture medium (DMEM-high glucose con-taining 10% fetal bovine serum (FBS) and 1% penicillin-
FIG. 1.
Schematics of proposed mechanism of tissue failure.
Healthy cells (green) secrete factors (grey) that enable or enhancethe function of other cells. Advective flow transport cooperative fac-tors downstream. Once a cell fails/dies (red), it can no longer supportcells that are dependent on it. A cascade of failure then propagatesover time in the direction of flow. streptomycin). The photoinitiator Irgacure D-2959 (Sigma-Aldrich) (final concentration: 0.1%, w/v) was added.
Cell culture and seeding.
Cardiac fibroblasts were isolatedfrom hearts of 2-day-old Sprague-Dawley rats (Charles RiverLaboratories), according to the IACUC guidelines with theapproval of the University of Notre Dame, which has an ap-proved Assurance of Compliance on file with the National In-stitutes of Health, Office of Laboratory Animal Welfare. Ratswere sacrificed via CO treatment, and the hearts were imme-diately collected, minced and incubated in trypsin (Life Tech-nologies) at 4 °C for 16h with gentle agitation as describedpreviously [13]. After further digestion with collagenase typeII (Worthington-Biochem) at 37°C, the tissues were strainedthrough 40 µ m filters and the cells in filtrate were incubatedat 37°C for 2h. This brief incubation allowed exclusively fi-broblasts to attach on the plate. After removal of unattachedcells, Fibroblasts were incubated in culture media with me-dia change every 3 days until passage 4 (P4). When the cellswere at approximately 80% confluency, they were detachedfrom flasks using 0.25% trypsin-EDTA (Life Technologies),counted, and reconstituted in culture medium.The cell suspension was mixed with the PEG:PEG-RGDsolution at 1:1 volume ratio, and loaded into the microfluidicdevice such that the final cell density was × cells/mL,PEG:PEG-RGD concentration was 10% (w/v), and photoini-tiator concentration was 0.05% (w/v). The polymer was ex-posed to UV at 365nm wavelength and 6.9mW/cm intensityfor 60s to crosslink it. The device was connected to a pumploaded with a syringe full of culture medium, and perfused at aflow rate of µ L/h for 30min to wash away the photoinitia-tor and the uncrosslinked polymer remaining in the gel. Thebright field images of the gel were taken using a microscope(Zeiss, Hamamatsu ORCA flash 4.0), and the initial total cellnumber was determined using Fiji software (NIH).
FIG. 2.
Time-lapse imaging of the cardiac fibroblast-laden PEG:PEG-RGD gels in the microchannel showing cell viability in real-time.(a)
Time-lapse fluorescence microscopy images showing the dead cells in the gel every 2 h for 16h. Red: ethidium homodimer-1. (b, c)
Linegraphs showing the change in cell viability along the microchannel over time. (b)
Dead cell number, and (c) cell viability. Flow rate: 20 µ L/min.
Cell viability measurements
Cell viability was determined using the Live/Dead cell via-bility assay (Thermo Fisher Scientific) either daily (at Days 0,1, and 2) or in real-time (every 2h for 16h). For the daily anal-ysis, the microfluidic device was disconnected at each timepoint, perfused with PBS containing ethidium homodimer-1(EthD-1) (4 µ M) for 30min, and imaged under bright field (to-tal cells) or fluorescence (dead cells, red) modes with a flu-orescence microscope (Zeiss, Hamamatsu ORCA flash 4.0).The device was reconnected to the syringe pump full of cul-ture medium, and perfused at µ L/h flow rate until the nexttime points (Day 1 and Day 2).To monitor cell viability in real-time, the device was per-fused with culture medium containing EthD-1 (4 µ M) at twoflow rates (20 and µ L/h), and imaged under bright fieldat time 0 (total cells), and then under fluorescence with time-lapse imaging at every 2h for 16h (dead cells). The cell num-bers in gel areas of 0.5mm × EXPERIMENTAL RESULTS
To investigate how failure propagates through a tissue un-der flow conditions, a microfluidic channel filled with cardiac fibroblast-laden hydrogel was perfused at various flow ratesand imaged in real-time every 2h for 16h. At 20 µ L/min flowrate, we observed a gradual increase in the dead cell numberover time in the inlet, while dead cell number remained rela-tively stable over time at the outlet (Fig. 2a,b). The differencebetween the dead cell number at t =
2h and at t =
16h gradu-ally decreased towards the outlet. Thus, after 16h of perfusion,cell viability was significantly higher at the outlet portion rel-ative to the inlet portion (Fig. 2c).To see the effect over a long time, in another experiment,the engineered tissue was perfused at 20 µ L/h flow rate for 2days, and the dead cells were imaged and counted at Days0, 1 and 2 (Supplementary Fig. 1). Upon perfusion, morecells died at the inlet and cell viability showed an increasingtrend along the gel towards the outlet. Viability at the inletwas significantly lower than that at the outlet both for Day 1( p < . ) and for Day 2 ( p < . ).Considering that the flow rate did not change along the mi-crochannel, the shear stress exerted on the cells by the flowitself should be the same. Thus, we concluded that the gradualincrease in cell viability towards the outlet was not because ofthe shear force, but rather due to accumulation of the cooper-ative factors.The difference between the relative cell viability (viabilityat a given time relative to that at t=2h) at the inlet and outletchanged dramatically when the flow rate was altered (Fig. 3).At 20 µ L/min flow rate, relative cell viability decreased by 2%at the outlet and 14% at the inlet over a 16h period, with agradual decrease at both edges (Fig. 3a). The viability differ-ence between the inlet and outlet reached 12% at t = µ L/min flow rate, viability at the outlet decreased by 58%
FIG. 3.
Experimental data with fitted curves.
Fitted parametervalues are k = 1 . , α = 0 . , β inlet = 5 . , β outlet = 7 . for v = µ L/h, and β inlet = 1 . × − , β outlet = 0 . for v = µ L/h.We indeed see β inlet ¡ β outlet for each flow rate, signifying a larger por-tion of cooperative factors in the vicinity of cells at the outlet. Inletcell populations were given by averaging over measurements takenbetween 1.25 and 2.5 mm from the inlet for v = µ L/h and 3.75and 5 mm from the inlet for v = µ L/h. Outlet cell populationswere given by averaging over measurements taken between 5 and6.25 mm from the inlet for v = µ L/h and 8.75 and 10 mm fromthe inlet for v = µ L/h. Error bars correspond to one standarddeviation from the mean. over a 16h period, while that at the inlet decreased by 93%(Fig. 3b). Although cell viability was initially higher at theinlet (69%) than the outlet (38%), it decreased dramaticallyover time and became higher at the outlet (16%) than the inlet(11%).Therefore, at the low flow rate, cell viability was the high-est, but the difference between the viability at the inlet andoutlet was the smallest. With the increasing flow rate, viabil-ity decreased remarkably especially at the inlet but also at theoutlet. These observatons suggest that, at lower flow rates,cooperative factors were not completely cleared from the inletand could bind to their receptors on cell surfaces, improvingcell viability outcomes. When the flow rate was increased,more cooperative factors were washed out from the inlet, giv-ing them less time to bind to receptors on the cells and dra-matically decreasing cell viability. As cells at the inlet dieout, they are no longer able to produce cooperative factors toaid the survival of cells at the outlet, leading to a propagationof failure and a reduced viability of cells at the outlet.
ANALYTICAL RESULTSModel
To characterize the dynamics of flow mediated failure prop-agation in a more general setting, and to make predictions be-yond the experiments described above, we develop an ana-lytical model the inputs of which we obtain from experimentsdescribed above. Our model assumptions, qualitatively stated,are as follows. (1) The cells do not proliferate and migrate, butare under stress and die. We denote the cell density at position x along the channel at time t , with n ( x, t ) . (2) Cells secretesome cooperative factor(s) that will diffuse, flow and decaywithin the circulating fluid, and help other cells survive/func- tion. We assume that the diffusion and decay parameters forthese factors are similar, and denote their concentration col-lectively by a single quantity Φ( x, t ) . These assumptions canbe quantified as ˙ n = − α φ k φ k + Φ k n (1) ˙Φ = d ∇ Φ − v · ∇ Φ − γ Φ +
An. (2)where, d is the diffusion constant for the cooperative factors, v is the flow velocity, A is the rate at which cells secrete co-operative factors, k describes the steepness of the responseof cells to cooperative factors, and the constant φ quantifiesthe “required amount” of factors for a cell to function nor-mally. The functional form of the right hand side of eqn.1 isthe Hill function, which accurately describes a cell’s responseto many different kinds of molecular agents. It is motivatedby Michaelis-Menten type reaction kinetics, fits experimen-tal findings, and is ubiquitously used in population dynamicsmodels. α is the cell death rate when there are no cooperativefactors. In a variable environment α could be time dependent,however here we take it to be constant.We study the phenomena of failure propagation by first ex-tracting the relevant length scales. From eqn.(2), we can solvethe Green’s function for a point source (Appendix II). Fromthe Green’s function, we can then extract left (cid:96) L and right (cid:96) R length scales given as, (cid:96) L ( v ) = (cid:104) λ ( v ) + v d (cid:105) − , (cid:96) R ( v ) = (cid:104) λ ( v ) − v d (cid:105) − , (3)where we have defined λ ( v ) = (cid:112) γ/d + v / d and flow istaken to be from left to right ( v ≥ . These length scales cor-respond to the characteristic advection-diffusion-decay lengthof the cooperative factors given by a point source.As the flow velocity v increases, the left length goes to zeroand the right length increases roughly linearly with respect to v . This introduces a bias in the concentration of cooperativefactors, and we expect to see a propagation of cell death inthe direction of flow. As cells upstream die, the cooperativefactors upstream also diminish, causing a propagation of deathmoving downstream, in the direction of flow. Experimental fits
To fit our model to experimental data, we first make somesimplifying assumptions, in order to reduce our model to asolvable system. We assume the local cooperative factor con-centration Φ( x, t ) will be proportional to the local cell density n ( x, t ) . We denote this proportionality constant β , so that Φ( x, t ) = βn ( x, t ) . In the center of the microchannel, where n is roughly uniform, this constant would be given by tak-ing the steady state cooperative factor concentration, giving Φ =
An/γ , so β = A/γ . By the inlet, this constant willbe lower in value since the area of cells contributing to thecooperative factors will be reduced.The constant β will in general depend on the length scalescontributing to the local cooperative factor concentration and FIG. 4.
Simulation results for various regimes.
If the initial cellconcentration is sufficiently large ( n > γφ /A ), there will be apropagation of death due to flow. Here we numerically solve and plotcell concentration normalized by the initial concentration at varioustimes normalized by the death rate α . (a) For low hill constant ( k =1 ), we see a permanent gradient formed before the bulk goes belowthe critical value, n < γφ /A , and cells begin to die uniformly. (b) As we increase the hill constant ( k = 2 ), the bulk cell densityattenuates slower and a more pronounced wave of failure begins todevelop. (c) In the case of a strong response function ( k = 4 ), wesee a more pronounced wave of propagating failure. Here, the bulkconcentration remains roughly constant and death propagates fromthe inlet to outlet, along direction of flow, at roughly a constant deathvelocity. thus will vary between the inlet and outlet regions. Specifi-cally, we expect the value of β inlet to scale as A(cid:96) L / ( (cid:96) L + (cid:96) R ) γ and β outlet ∼ A(cid:96) R / ( (cid:96) L + (cid:96) R ) γ . Therefore, for a large flowrate, we expect the value of β inlet to be less than β outlet sincemore cells are able to contribute to the cooperative factor con-centration at the outlet. With this assumption, we can writedown an effective equation describing the growth of the cellsover time as, ˙ n = − αn βn/n ) k , where n is the initial cell concentration. We can then solvethis exactly to get, n ( t ) = ( n /β ) (cid:104) W (cid:16) e β k − αkt β k (cid:17)(cid:105) /k where W ( z ) is the Lambert W function, defined as the prin-ciple solution for w in z = we w and can be computed toarbitrary numerical precision.We then fit our model to experimental data in Fig.3. Sincewe expect the growth kinetics to be the same in both experi-ments, and the only difference being flow, we constrain the fitparameters for α and k to be the same for all fits and allowdifferent values of β for inlet and outlet regions and for eachflow rate v = µ L/h and v = 120 µ L/h. For hill and decayconstants we then get k = 1 . , α = 0 . − . For flowrate v = µ L/h, we get β inlet = 5 . , β outlet = 7 . . For v = 120 µ L/h, we get β inlet = 1 . × − , β outlet = 0 . . We see that indeed β inlet < β outlet as expected for each flowrate. We also see the values of β are much smaller for thelarger flow rate, with the inlet value of β being vanishinglysmall at the flow rate v = 120 µ L/h. This is because with withlarge flow rate, the cooperative factors are pushed much fur-ther downstream and no longer help the cells at the inlet andmay also diminish the cooperative factor concentration at theoutlet.We next study the system of equations (1-2) and derive con-ditions for a propagation of failure, as well as the velocity andacceleration of failure propagation.
Model regimes and failure propagation
From numerical simulations and dimensional analysis, wefind the condition for failure propagation to occur is that theinitial density of cells n must be sufficiently large, such that An /γ (cid:29) φ . This is since the largest the cooperative factorconcentration can be is given by An /γ and this value mustbe above the threshold concentration φ for cell density to notdecay exponentially. If the cell concentration drops below thiscritical value, there will no longer be a pronounced propaga-tion of failure. Instead, the entire population of cells will dieroughly uniformly at an exponential rate of α .For cases where there is a propagation of failure, where An /γ > φ , we plot the numerical solution to our systemin Fig. 4, for values of the hill constant k = 1 , , . We findthat as we increase the value of the hill constant, we get a morepronounced wave. The population dynamics is strongly deter-mined by the initial density n and the form of the responsegiven by the hill constant k . In the case where k is low, thebulk where Φ( x ) > φ will attenuate quicker. We find thisattenuation of the bulk leads to an “acceleration” of death. Wealso study this attenuation and the corresponding accelerationof death propagation.We now determine the velocity at which the death of cellswill propagate. Our derivations for failure propagation veloc-ity, depth, and acceleration are given in detail in Appendix III.We illustrate here the procedure used in our derivations.To simplify our analysis, we used a boxcar approxima-tion to the Green’s function for the cooperative factors whereleft and right lengths are given by equations (3), and assumechemicals quickly reach steady state. The total area of theGreen’s function corresponds to the secretion rate per cell den-sity at steady state. We therefore have G ( x ) = ( A/γ )Θ( x + (cid:96) L )Θ( (cid:96) R − x ) , where Θ( x ) is the Heaviside step function.We then convolve this boxcar Green’s function with a semi-infinite initial cell concentration, n ( x,
0) = n Θ( x ) . Thisgives a cooperative factor concentration profile that increaseslinearly up until x = (cid:96) R , after which Φ( x ) saturates to a con-stant given by Φ( x > (cid:96) R ) = An /γ . Assuming the initial celldensity is sufficiently large such that An /γ > φ , we canfind the position ∆ , such that Φ(∆) = φ . To the left of thisvalue, cells are expected to decay exponentially at a rate α .The cooperative factor concentration will then update to thisnew concentration of cells and we can reiterate this to get thenext decay of cells. The death of cells will therefore continuepropagate by an amount ∆ at a rate α , giving a first approx-imation to the failure propagation velocity as v d = α ∆ . Wefind this gives good agreement for large k , where the cell deathrate behaves closer to a step function, and for short times.We then further improve these calculations by performinga second iteration with an updated approximate cell concen-tration, taking into account the region where φ < Φ( x )
Analytical theory and numerical simulations for failurepropagation and bulk death versus time. (a)
Bulk cell concentra-tion normalized by initial cell concentration n b /n , versus time. Forlower values of the hill constant k , the cell concentration in the bulkdecays quicker. At a critical time t c , the bulk density falls below thecritical value n b < γφ /A and dies out exponentially. (b) Propaga-tion of failure versus time for hill constants k = 1 , , , . For lowerhill constants, the failure depth rapidly increases up to a critical time t c , where the bulk collapses before the propagating wave reaches theend of the domain. For larger hill constants, the velocity remainsroughly constant and reaches the end of the bulk before time t c . (c) Failure propagation velocity vs flow velocity. The failure propaga-tion velocity v d increases roughly linearly with the flow velocity v . (d) Failure propagation acceleration vs hill constant. As the hill con-stant k increases, the failure propagation becomes more pronouncedand moves at a more constant speed. The acceleration then decreasesroughly exponentially as k increases. Solid lines are obtained fromnumerical simulations throughout and dotted lines are obtained fromanalytical theory. DISCUSSION
Here we studied how failure propagates in a system whereinterdependence is mediated by flow. These results emphasizethe importance of intercellular processes on aging.We performed experiments with synthetic tissues filled withcardiac fibriblast-laden PEG hydrogels in a microchannel andobserved that flow can help increase the lifespan of cellsdownstream of the flow (Fig.1, 2). We explained this observa-tion with cooperative factors, which were carried by the flowtowards the outlet. Cooperative factors are known to promotecell survival [13–20], cardioprotection [15, 16], and angiogen-esis [16, 17], and cells failing to receive the necessary factorsfrom the neighboring cells go through apoptosis [25]. Moti-vated by these results, we developed an analytical model todescribe the death of cells that communicate via diffusive co-operative factors in a flowing environment.Fitting this model to experiment, we saw indeed that theproportion of cooperative factors downstream of the flow werelarger than those upstream. This then leads to a faster death ofcells by the inlet and longer lifespan for cells by the outlet.Next, we investigated further the consequences of our an-alytical model. We found analytical conditions for a “wave”of failure propagation to occur in the direction of flow. Ascells die out upstream, this failure will propagate and increasethe mortality rate of cells downstream. The conditions for thispropagation to occur were found to be a sufficient density ofcells and a non-zero flow rate. Once cell density decreases be-low a critical threshold, the cells will die uniformly at roughlythe same exponential rate.In the case of sufficiently large cell density, we see the formof failure propagation is heavily dependent on the hill constant k . For smaller values of k , cells that are further downstreamdie faster compared to large k , leading to an “acceleration”of failure. For large constants k , the wave front is more pro-nounced and continues at more of a constant velocity.Through a dimensional analysis and simplifying assump-tions, we were able to derive analytical formulas describingthe propagation depth, velocity, and acceleration of failure, aswell as the attenuation of cells beyond the wave front.Failure propagation has already been studied previouslywhen cells are coupled through diffusion [13, 14]. Here, weanalyzed the failure propagation under flow, and developed ananalytical model that describes the depth, velocity, and accel-eration of failure propagation. This is physiologically more relevant, since cooperative factors are carried via flow in thebody. For example, paracrine factors are transported via inter-stitial flow [26, 27].Fluid flow is known to mediate communications within mi-crobial communities and influence death patterns [28]. Thus,while our model was motivated by experiments on mam-malian tissues, we might expect similar results to also holdfor eukaryotic colonies and bacterial biofilms where diffusiveand advective forces are responsible for the communication ofinterdependent members of the population. AUTHOR CONTRIBUTIONS
GU, GB, PZ and DCV formulated the problem. GB andPZ designed and carried out the experiments. GU and DCVdesigned and carried out the theory. GU, GB, DCV wrote thepaper.
ACKNOWLEDGMENTS
This study was funded by National Science Foundationgrant CBET-1805157 [1] S. I. Rattan, Free radical research , 1230 (2006).[2] E. H. Blackburn, Nature , 53 (2000).[3] J. M. Murabito, R. Yuan, and K. L. Lunetta, Journals of Geron-tology Series A: Biomedical Sciences and Medical Sciences ,470 (2012).[4] J. W. Shay and W. E. Wright, Nature reviews Molecular cellbiology , 72 (2000).[5] L. A. Gavrilov and N. S. Gavrilova, Journal of theoretical Biol-ogy , 527 (2001).[6] D. C. Vural, G. Morrison, and L. Mahadevan, Physical ReviewE , 022811 (2014).[7] J. W. Vaupel, J. R. Carey, K. Christensen, T. E. Johnson, A. I.Yashin, N. V. Holm, I. A. Iachine, V. Kannisto, A. A. Khazaeli,P. Liedo, et al. , Science , 855 (1998).[8] G. Caughley, Ecology , 906 (1966).[9] N. Stroustrup, W. E. Anthony, Z. M. Nash, V. Gowda,A. Gomez, I. F. L´opez-Moyado, J. Apfeld, and W. Fontana,Nature , 103 (2016).[10] S. G. Farrell, A. B. Mitnitski, K. Rockwood, and A. D. Ruten-berg, Physical Review E , 052409 (2016).[11] A. Mitnitski, A. Rutenberg, S. Farrell, and K. Rockwood,Biogerontology , 433 (2017).[12] J. J. Boonekamp, M. Briga, and S. Verhulst, Experimentalgerontology , 95 (2015).[13] A. Acun, D. C. Vural, and P. Zorlutuna, Scientific reports , 1(2017).[14] D. Suma, A. Acun, P. Zorlutuna, and D. C. Vural, Royal Soci-ety open science , 171395 (2018).[15] L. Zhao, X. Liu, Y. Zhang, X. Liang, Y. Ding, Y. Xu, Z. Fang,and F. Zhang, Experimental cell research , 30 (2016).[16] H. Li, S. Zuo, Z. He, Y. Yang, Z. Pasha, Y. Wang, and M. Xu,American Journal of Physiology-Heart and Circulatory Physi-ology , H1772 (2010). [17] S. Jiang, H. K. Haider, N. M. Idris, A. Salim, and M. Ashraf,Circulation research , 776 (2006).[18] C.-C. Lin and K. S. Anseth, Proceedings of the NationalAcademy of Sciences , 6380 (2011).[19] M. Mil´an, L. P´erez, and S. M. Cohen, Developmental cell ,797 (2002).[20] F. Boissard, M. Tosolini, L. Ligat, A. Quillet-Mary, F. Lopez,J.-J. Fourni´e, L. Ysebaert, and M. Poupot, Oncotarget , 52225(2017).[21] C. A. Davison, S. M. Durbin, M. R. Thau, V. R. Zellmer,S. E. Chapman, J. Diener, C. Wathen, W. M. Leevy, and Z. T.Schafer, Cancer research , 3704 (2013).[22] S. Fonseca, M. Reis, V. Coelho, L. Nogueira, S. Monteiro,E. Mairena, F. Bacal, E. Bocchi, L. Guilherme, X. Zheng, et al. ,Scandinavian journal of immunology , 362 (2007).[23] I. Tamm and T. Kikuchi, Journal of cellular physiology ,494 (1990).[24] O. R. Jones, A. Scheuerlein, R. Salguero-G´omez, C. G. Ca-marda, R. Schaible, B. B. Casper, J. P. Dahlgren, J. Ehrl´en,M. B. Garc´ıa, E. S. Menges, et al. , Nature , 169 (2014).[25] R. S. Petralia, M. P. Mattson, and P. J. Yao, Ageing researchreviews , 66 (2014).[26] W. Yao, Y. Li, and G. Ding, Evidence-Based Complementaryand Alternative Medicine (2012).[27] M. A. Swartz and A. W. Lund, Nature Reviews Cancer , 210(2012).[28] T. Rossy, C. D. Nadell, and A. Persat, Nature communications , 1 (2019). APPENDIXI. ADDITIONAL EXPERIMENTAL RESULTS
To see the effect over long time, the engineered tissue wasperfused at 20 µ L/h flow rate for 2 days, and the dead cellswere imaged and counted at Days 0, 1 and 2 (SupplementaryFig. 1a(i)). At Day 0, bright field images were also takento account for the initial total cell number. Cell viability atDay 0 was similar throughout the gel, except the viability atthe inlet (point 1.25 mm) being significantly higher than theother points ( p < . when compared to outlet (point 7.5mm)) (Supplementary Fig. 1a(ii, iii)). Upon perfusion, morecells died at the inlet and cell viability showed an increasingtrend along the gel towards the outlet. Viability at the inletwas significantly lower than that at the outlet both for Day 1( p < . ) and Day 2 ( p < . ).In another experiment, we kept all the parameters the same,but reversed the flow direction after Day 1 (SupplementaryFig. 1b(i)). At Day 0, cell viability was the same all along thegel. Upon perfusion for 1 day, more cells died at the inlet (leftside of the gel) than the outlet (right side), and cell viabilityfollowed an increasing trend towards the outlet (Supplemen-tary Fig. 1b(ii)). After day 1, flow direction was changed;perfusion was applied from the right side of the gel. After per-fusion for another day, but in the reverse direction, we againobserved an increasing cell viability towards the outlet (leftside of the gel). While, at Day 1, cell viability was signifi-cantly higher at the right side of the gel (outlet) than the left( p < . ), it was higher at the left side than the right atDay 2 ( p < . ) (Fig. 1a(iii)). II. CHEMICAL GREEN’S FUNCTION
We assume the simplified case of a constant flow profile, v = v ˆ z . Assuming a constant flow profile, we can use a shiftof coordinates to simplify the chemical dynamics. If we let x d be the coordinate in the direction of flow, and let z = x d − vt ,the chemical equation reduces to ∂ Φ ∂t = D ∇ Φ − γ Φ + An where now, ∇ = d − (cid:80) i =1 ∂ x i ˆ x i + ∂ z ˆ z .Now for an arbitrary source function S ( x , t ) = An ( x , t ) ,the chemical profile will be given by, Φ( x , t ) = G ( x , t ) ∗ S ( x , t ) , (1)where ∗ is the convolution operator. The Green’s function G ( x , t ) is the solution to (cid:0) ∂ t − D ∇ + γ (cid:1) G ( x , t ) = δ ( x ) δ ( t ) Taking a Fourier transform with respect to space ( x ) , we get (cid:0) ∂ t + k D ∇ + γ (cid:1) ˜ G ( k , t ) = δ ( t ) . The solution to this is given by, ˜ G ( k , t ) = Θ( t ) exp (cid:2) − ( k D + γ ) t (cid:3) . We recognize this as the Green’s function for the diffusionoperator times a decay factor of exp[ − γt ] . We then get, forthe Green’s function in d dimensions, G ( x , t, v ) = Θ( t ) (cid:18) πDt (cid:19) d/ e − r / Dt e − γt (2)where r = d − (cid:80) i =1 x i + ( x d − vt ) .If we now have a stationary point source, we can get thesteady state Green’s function by convolving with a source, S ( x, t ) = δ ( x ) constant in time.In the one dimensional case, we get, Φ( x, v ) = ∞ (cid:90) −∞ ∞ (cid:90) −∞ G ( x (cid:48) − x, t (cid:48) ) δ ( x (cid:48) ) d x (cid:48) d t (cid:48) = exp (cid:104) x (cid:16) v d − sign( x ) (cid:112) γ/d + v / d (cid:17)(cid:105) d (cid:112) γ/d + v / d ≡ exp (cid:104) x (cid:16) v d − sign( x ) λ (cid:17)(cid:105) / (2 dλ ) where we have defined λ ( v ) = (cid:112) γ/d + v / d . From here,we can get diffusion-advection-decay lengths given by (cid:96) L ( v ) = (cid:104) λ ( v ) + v d (cid:105) − (cid:96) R ( v ) = (cid:104) λ ( v ) − v d (cid:105) − , (3)assuming the flow is from left to right ( v ≥ . III. FAILURE PROGAGATION VELOCITY DERIVATION
We begin by assuming that the cells respond to the chemicalconcentration “sharply”. This is the case as k → ∞ , wherethe response function becomes a step function. We then laterrelax this assumption.To simplify our analysis, we used a boxcar approximationto the Green’s function for the cooperative factors where leftand right lengths are given by eqn.3, and assume chemicalsquickly reach steady state. The total area of the Green’s func-tion corresponds to the secretion rate per cell density at steadystate. We therefore have G ( x ) = ( A/γ )Θ( x + (cid:96) L )Θ( (cid:96) R − x ) , (4)where Θ( x ) is the Heaviside step function.To get a first approximation to the cooperative factor con-centration, we convolve the boxcar Green’s function (equa-tion 4) with an initial semi-infinite concentration of cells, n ( x,
0) = n Θ( x ) . This gives us a linearly increasing chem-ical profile, with slope H/W , that saturates to a maximum
Supplementary Figure 1.
The change in viability of cardiac fibroblasts along the microchannel upon perfusion.
Cell viability in thePEG:PEG-RGD gel perfused for (a) two days from left to right or (b) one day from left to right followed by one day from right to left. (i)Fluorescence microscopy images showing the dead cells in the gels at Day 0, Day 1 and Day 2. Red: ethidium homodimer-1. (ii) Line graphsshowing cell viability along the gels. (iii) Bar graphs showing cell viability at the inlet and outlet portions of the gels. Flow rate: 20 µ L/min.Asterisks indicate statistical significance p < . value of H after a length of (cid:96) R into the microchannel, where H = An /γ and W = (cid:96) L + (cid:96) R . Specifically, we get Φ( x ) = , x < − (cid:96) LHW ( x + (cid:96) L ) , − (cid:96) L < x < (cid:96) R H, x > (cid:96) R (5)Now, if the response to the chemicals is a step function (asis the case when k → ∞ ), and if H (cid:29) φ , cells where Φ( x ) is below φ will die at a rate α and cells where Φ( x ) is above φ will survive. The cell density will then shift to the right byan amount ∆ given by setting Φ(∆) = φ or, ∆ = 12 Aγn (cid:104) An v + (2 γφ − An ) (cid:112) dγ + v (cid:105) . (6)The cooperative factor concentration will then update withthis new concentration of cells and we can then reiterate thisto get the next death of cells. The death of cells will thencontinue to propagate by an amount ∆ at a rate α . The deathpropagation velocity is therefore given as v d = α ∆ . We findthis gives good agreement for large k and short times, but failsto capture the strong time dependence for small k .We then improve on our analytical results by consideringthe death of the bulk (in the region x > (cid:96) R ). We replace theinitial cell concentration n in equation (6) by a time varyingfunction n b ( t ) for the bulk cell concentration. To get n b ( t ) ,we approximate the chemical concentration in the bulk as thesteady state, non-spatial concentration, Φ ∗ = An b ( t ) /γ . Wethen expand the hill form about infinity (high concentration limit), to get /φ ) k + 1 ≈ (cid:18) φ Φ (cid:19) k − (cid:18) φ Φ (cid:19) k + (cid:18) φ Φ (cid:19) k − . . . Taking the first order term and plugging in Φ ∗ , we get for theattenuation of the bulk cell density, n b ( t ) = (cid:34) n k − α (cid:18) γφ A (cid:19) k kt (cid:35) /k . (7)Plugging this into the above expression for v d = α ∆ , we get, v d = αv γ + α (cid:112) dγ + v γ (cid:34)(cid:18) An γφ (cid:19) k − αkt (cid:35) − k − . This expression for the velocity better captures the accelera-tion of velocity over time, but grows in error over time forsmall k . This is because we fail to consider the region where Φ > φ but not yet fully saturated. This gives another re-gion where n ( x, t ) does not die out exponentially, but dies outfaster than n b ( t ) and will lead to an increase in the velocity v d .Therefore, for a better approximation to the velocity, wealso take into account the region where the cooperative factorconcentration is larger than the threshold φ , but not yet satu-rated. The cell concentration in this region (denoted n δ ( x, t ) )will decay faster than the cell concentration in the saturatedregion, n b ( t ) .0 Supplementary Figure 2.
Schematic of theoretical methods. (a)
First iteration. The initial semi-infinite concentration of cells resultsin a piecewise linear concentration of the cooperative factors. Theposition where the cooperative factors equal the threshold value φ isgiven by ∆ . Below this value, cells will die out exponentially. A firstapproximation to the failure propagation velocity can then be givenby v d = α ∆ . (b) Second iteration. To improve our calculation of thefailure propagation velocity, we account for the fact that the regionbetween ∆ and (cid:96) R in the first iteration dies out faster than the regionafter (cid:96) R . We approximate this intermediate region by a flat spatiallyconstant value that is less than the cell concentration after (cid:96) R . Wethen solve again for when the cooperative factor concentration equalsthe threshold value φ and get a larger value for ∆ . We then approximate the cell concentration n ( x, t ) as n δ ( t ) in the region where φ < Φ( x ) < An b /γ and n b ( t ) where Φ( x ) > An b /γ , as shown in Fig. 2.Taking the convolution of this updated cell concentrationwith our boxcar Green’s function, we get another piecewisefunction for Φ( x ) , with 5 regions. Our second iteration of Φ( x ) then gives, Φ = x < − (cid:96) LH δ W ( x + (cid:96) L ) − (cid:96) L < x < − (cid:96) L + δ H ∞ W ( x + (cid:96) L ) − δ ( H ∞ − H δ ) W − (cid:96) L + δ < x < (cid:96) R H δ + H ∞ − H δ W ( x + (cid:96) L − δ ) (cid:96) R < x < (cid:96) R + δH ∞ x > (cid:96) R + δ where H δ = An δ /γ and H ∞ = An b /γ . We now solve againfor Φ(∆) = φ to get a new expression for ∆ . We assume δ issmall, and use the center condition to solve for the thresholdcrossing. This assumption works for large velocities as ∆ → (cid:96) R and δ = (cid:96) R − ∆ → . Solving for the threshold crossingthen gives, ∆ = ( W φ − H ∞ ( (cid:96) L − (cid:96) R ) − H δ (cid:96) R ) / (2 H ∞ − H δ ) ,which in terms of original paramters gives, ∆ = Av (2 n b − n δ ) + (cid:112) v + 4 dγ (2 γφ − An δ )2 Aγ (2 n − n δ ) , (8)Note, that this reduces to our original formula if we set n δ = n b as expected.Now, to get the velocity propagation, we need an expres-sion for n δ . We could in principle obtain this from solving for n ( x, t ) with a linear approximation for Φ( x ) . We instead usea simplified expression from the following arguments. As thehill constant k → ∞ , the cell concentration will remain con-stant for any value of Φ above the threshold φ , and the two regions should remain the same, n δ ( t ) = n b ( t ) . For k = 1 weapproximate n δ ≈ n b . We therefore approximate n δ ( x, t ) as n δ ( t ) = (cid:34) − (cid:18) (cid:19) k (cid:35) n b ( t ) . (9)Now taking the velocity to be v = α ∆ and substituting equa-tions (7) and (9) for n b and n δ , we get, v d = αv γ + α (cid:112) v + 4 dγ k ) × − k γ + 2 k φ A (cid:34) n k − αkt (cid:18) γφ A (cid:19) k (cid:35) − k (10)We find this gives good agreement with numerical simulationsof the model.We can also integrate this to get the failure penetrationdepth over time, x d ( t ) = (cid:82) t v d ( t (cid:48) ) d t (cid:48) , giving, x d = (cid:112) v + 4 dγ k ) Aγ (cid:40) At (cid:32) (1 + 2 k ) v (cid:112) v + 4 dγ − (2 k − (cid:33) + 2 k γφ α (1 − k ) (cid:18) γφ A (cid:19) − k (cid:32) n k − αkt (cid:18) γφ A (cid:19) k (cid:33) k − k − n k − for k (cid:54) = 1 , and x d = vt γ − (cid:112) v + 4 dγ αγ (cid:26) αt + 4 log (cid:20) − αt (cid:18) γφ An (cid:19)(cid:21)(cid:27) for k = 1 . Also, taking the derivative of equation (10) givesthe acceleration, a d = 2 k αφ (cid:112) v + 4 dγ (cid:16) γφ A (cid:17) k (cid:20) n k − αkt (cid:16) γφ A (cid:17) k (cid:21) − k +1 k (1 + 2 k ))