Interdependence theory of tissue failure: bulk and boundary effects.

The mortality rate of many complex multicellular organisms increases with age, which suggests that net ageing damage is accumulative, despite remodelling processes. But how exactly do these little mishaps in the cellular level accumulate and spread to become a systemic catastrophe? To address this question we present experiments with synthetic tissues, an analytical model consistent with experiments, and a number of implications that follow the analytical model. Our theoretical framework describes how shape, curvature and density influences the propagation of failure in a tissue subjected to oxidative damage. We propose that ageing is an emergent property governed by interaction between cells, and that intercellular processes play a role that is at least as important as intracellular ones.

Introduction

As an organism ages, its cells shrink or enlarge, increase their lipid and pigment content, and lose their functionality and ability to proliferate. Mechanical theories of ageing typically focus on the biomolecular mechanisms governing metabolism, cell damage and repair, such as oxidative stress, shortening telomeres and various genetic factors [1-5].

However, ageing is a complex phenomenon that spans multiple time and length scales. Organisms do not ultimately die because they run out of cells, but because cellular damage manifests as larger-scale problems in tissues and organs, through a cascade of interactions. It remains unclear how failures dynamically propagate and accumulate to lead to frailty, ageing-related diseases, and ultimately, death [6-10].

The catastrophic end state, universality of survival statistics, and the irreversible nature of failure can be understood in terms of the dynamics of a network of interdependent components [6]. According to this interdependence network theory of ageing, if a component fails, then any other component that crucially depends on it will also fail. As a result, in well-connected interdependence networks, small failures will cascade into larger ones, and the probability of a system-wide catastrophe will monotonically increase with time, consistent with experimentally observed survival curves. The theory also predicts that the catastrophe is unavoidable even when repair rates far exceed the damage rate.

Some experimental work has been done, especially in cancerous tissues, to determine the effects of intercellular communication and stress [11-15]. These cooperative effects, our earlier experiments [9] and novel experiments reported here motivate us to view a tissue as a simple interdependence network where cells influence their neighbours’ response to damage via various diffusive factors.

Our aim here is to put forth a quantitative theory of tissue failure that is consistent with experimental data. In order to bridge microscopic (cellular) damage and macroscopic (organismic) catastrophes, we study the failure dynamics of an intermediate structure, the tissue. Specifically, we determine the relationship between damage propagation in functioning healthy tissues and the global and local properties, such as shape and density.

While the postulates of the quantitative model presented here are motivated by experiments on mammalian cells, we might expect similar results to also hold for eukaryotic colonies and bacterial biofilms where the survival of cells are linked to one other through a number of signaling and cooperative factors [16]. Thus, our results can be interpreted more generally, as the spatial dynamics of a cooperative population.

This paper is organized as follows: we first present experimental data on synthetic tissues that form the basis of our quantitative model. We then present analytical and computational results describing how failure propagates across a tissue.

Results

Motivating experiments

Rat fibroblasts were encapsulated within hydrogels of differing geometries, in particular, a triangular prism with 30°, 60° and 90° corners and a ‘flat’ cylinder. The distance between cells was measured to be quite uniformly distributed at approximately 40 μm, much larger than the size of a fibroblast, approximately 10 μm.

The cell-laden hydrogels were subjected to continuous oxidative stress through cell culture media containing hydrogen peroxide (0.2 mM). Throughout the duration of the culture, live and dead cells were stained and counted in order to measure the death rate at the edges, corners and bulk of the cell-laden hydrogels. The hydrogel material, which is arginine-glycine-aspartic acid (RGD) conjugated polyethylene glycol- 4-arm acrylate (PEG–RGD), was chosen to prohibit cells from migrating or proliferating [17,18]. This way we ensured that the observed effects were not due to a loss or gain of physical contact between cells, but due to intercellular cooperation via factors secreted by the cells that diffuse through the hydrogel medium. Furthermore, when the cell density is sufficiently low, the ageing effect disappears since the intercellular distance prohibits the interactions [9]. In this dilute limit the death rate increases, further supporting our hypothesis that ageing is driven by intercellular interactions rather than limitations in the externally provided resources.

We observed that cells located at smaller angles die faster than those at larger angles (figures 1 and 2). We also observed that the cells encapsulated in a triangular geometry died faster than those encapsulated in a disk geometry (figure 1a). In addition, cells started dying close to the edges and corners, and the dead layer propagated inwards. Specifically, when we plot the dead layer thickness as a function of time for the 30°, 60°, 90° and 180° cases, it was seen that the smaller angles experienced faster propagation than the larger ones (figure 1b).

Figure 1.

Dependence of population dynamics on geometry. (a) Images taken at time points 0, 8 and 16 h (top to bottom) after the start of stress treatment which show the collapse of each region of the geometry. Live cells (shown in green, first row) mark the initial conditions and dead cells (shown in red, remaining rows) show decay at the specified time points. A layer of dead cells forms and rapidly thickens until tips are completely consumed and the population collapse propagates inwards. Plotted in (b) is the penetration depth, i.e. the thickness of the dead cell layer, as a function of time for the 30°, 60°, 90° and 180° cases. We obtained the penetration depth by quantifying the dead cell number density at time 8 and 16 h (c), delineating the region separating the dead cell layer from the live cells (shown as thin white lines), and measuring the thickness of this layer from the tip. We used the dead cell number density as opposed to light intensity because edges appear brighter in the microscope images due to swelling of the edges, and not an increased dead cell population. The statistically significant differences have been marked with an asterisk, where a two tailed t-test with unequal variance gives p values of 1.6×10−4,0.012, and 0.006 for blue versus green, blue versus yellow, and yellow versus green curves, respectively, at the 16 h mark.

Figure 2.

Agreement of experiment (points) and theory (solid curve). The survival of cells acclimated in hydrogen peroxide for two days (a) and cells brought to Hayflick limit (b). Parameters used for the equation are β=3.8, α=1.5, k=3 (a) and β= 3.1, α=1.5, k=3 (b). For comparison purposes we have also included the behaviour of the null model where cells die independently of each other at a constant rate (dashed curve). Since the physiological requirements and biochemistry in both experiments should be similar, we have fixed k and α and obtained the fits by varying β using Mathematica's in-built χ2 test. The p values for the fits are 5×10−15 and 1×10−24. The moderate difference in the β value could be due to a difference between young and old cells in their cooperative factor secretion rates A, or their ability to receive/process cooperative factors Φ0.

These observations are unlikely to be caused by a diffusion limitation since the nutrient and oxygen diffusion penetration length is much larger than the height of the tissue, and thus oxygen and nutrients are in excess. Furthermore, the geometry of the samples are effectively two-dimensional planes, allowing every cell to receive equal amounts of external resources, which mostly diffuse from the top and bottom inward. Note that if the transport of nutrients and waste were the primary problem, we would expect to see the centre die first and propagate outwards, the opposite of what we observed.

The second, and perhaps more interesting reason that allows us to hypothesize that there is a cooperative interaction between cells is the following: if cells died only due to oxidation, then we would see a constant rate of death (e.g. 5% of the remaining live cells would die every day). However, we observe an accelerated rate of death (e.g. every day a larger percentage of the remaining live cells die). This suggests that lack of cells contributes an extra death rate in addition to oxidation and other single-cell-level causes of damage. This was also supported by our earlier experiments, where death rate was not time-dependent when cells were separated far enough that they could not interact [9].

It should be noted that the increased fluorescence around the edges of the hydrogels is due to swelling and induced edge curvature from the manufacturing process. Because of this, the fluorescence level was not used to count cells. Instead, only individual points were counted and remapped to density-based heat maps.

Analytical theory

To make sense of the geometry dependence of tissue lifetime and population dynamics we propose an analytical model, the assumptions of which, stated qualitatively, are as follows. (i) There exist damaging agents in the intercellular environment, and the cells counter this by secreting diffusive factors. Here we refer to the latter as ‘cooperative factors’ (CFs). The effect of CFs may be either direct or indirect, for example neutralizing oxygen species by directly binding to them, catalyzing or promoting intermediate byproducts that react with them, or participating in reactions or reaction cascades that activate higher order cell repair mechanisms. (ii) The probability of death of a cell increases with the relative abundance of the non-neutralized damaging agent, as given by first-order reaction kinetics. (iii) Cell lifetimes are much larger than the time it takes for the CF to diffuse throughout the sample, and the decay time of the molecule.

To be more precise, we assume that the concentration ϕ of CF at position r secreted by a single cell i located at r is approximately governed by the diffusion equation, where γ and D are the decay rate and diffusion constant of the CF in the extracellular matrix. The total concentration of CF experienced by cell i is then given by the sum of the secretions from all cells .

Secondly, given a local concentration Φ(r) near cell i, the probability of death per unit time is assumed to have a Michaelis–Menten–Hill form, where Φ0 is a constant that characterizes a threshold of cooperative factor concentration below which the cell’s survival is compromised, and k≥1 is the Hill constant, determined by the stoichiometry of the reactions in which ϕ is consumed; α is a proportionality constant that connects the probability of death to the concentration of cooperative factors or some unspecified molecule or structure that reacts with it.

We will now analytically solve the survival characteristics of the tissue in the bulk under certain reasonable approximations, and also simulate the system. The steady-state solution of (2.1), with appropriate boundary conditions, gives the influence of all surrounding cells on cell i: where n(,t) is the cell density at a given position and time. Every cell experiences a different amount of CF depending on the number of functional cells that surround it. In general, those near the boundary of the tissue have lesser neighbours. To calculate the overall loss of population, we calculate the cooperative factors received, averaged over position:

where N is the total number of cells. In the limit λL≪1, this expression is analogous to the problem of finding the electrostatic energy of a spherical charge distribution. In this limit, the result is well known, 〈ϕ(r)〉=(A/4π)6N/(5L). For an arbitrary λL≡d we obtain and for a uniform cell density N≡4πL3n(t)/3. As the cells die near the boundaries with higher likelihood, the uniformity assumption will break down. However, the agreement between discrete simulations and analytical theory indicates that the assumption of uniformity introduces a small overall error. We expect this error to be proportional to the surface area to volume ratio of the tissue.

The population of the cells is determined by their rate of death, dN/dt=−PN, which can be rewritten by substituting the approximation (2.5) into (2.2), This equation can be solved exactly: where N0 is the initial population of the tissue, and W(x) is the Lambert function, recursively defined by and can be expanded further by resubstituting W(x) on the right hand side, e.g. and so on.

In figure 2, we compare this theoretical result with experimental young and ‘old’ cell populations [9]. The former group was obtained by acclimating neonatal rat cells in oxidative media for two days, whereas the latter group was obtained by repeated cell division of the neonatal cells until they reached the Hayflick limit. The populations of both groups, encapsulated in a three-dimensional hyrogel, were measured every day. In both cases, we see that the qualitative shape of the survival curves are governed by the interdependence effect, although the old cells die quicker. The agreement between theory and experiment is acceptable.

In figure 3, we directly simulate the population dynamics of a tissue consisting of cells with random positions to compare it with our approximate analytical formula (2.8) and find good agreement.

Figure 3.

Dependence of population curves on diffusion length 1/λ (a), Hill coefficient k (b), and required cooperative factor threshold Φ0 (c). We compare simulations (dots) with analytical theory (solid line) for 1/λ=1,2,5,10,100 (black, blue, orange, green, red) units (a); and k=1,2,3,3.5 (black, blue, orange, green) (b), and Φ0=350,300,250,200 (green, orange, blue, black) (c). Remaining parameters are kept constant at N=2000, L=1, Φ0=300, k=3 and 1/λ=5 units.

To obtain the bulk lifespan τ of the tissue, we will first approximate the Lambert function by only keeping the leading term in its recursive definition, (the accuracy of this approximation is demonstrated in Fig. 4 with red dashed lines). This gives us,

where and C=αk/β (note that when β is large B∼1 and N(0)∼N0, as it should be). Then, from the approximate form (2.9), the lifetime τ can be estimated by setting n(τ)=0, Under the same approximation, the cell mortality rate (probability of death per unit time) is These approximate and exact analytical results are plotted in figures 3 and 4. The former shows the agreement between analytical theory and simulations.

Figure 4.

Population dynamics of a strongly interacting tissue model. In all panels the solid black curves are the exact analytical solutions, whereas the dashed curves are approximate solutions based on equation (2.9). In (a) and (b) we plot the fraction of surviving cells and cell death rate as a function of time. The system parameters are {α,k,N0}={1,2,1000}, and the solid curves, from thick to thin, correspond to β=2,4,6,8,10. In (c) and (d) we display the approximate formula for tissue lifetime as a function of β≈AN0/(Φ0λ2D), where A is the CF secretion rate, N0 the initial population, Φ0 the threshold CF level, 1/λ the CF decay length and D the diffusion constant of the CF in the medium. In (d) we plot failure propagation velocity from the surface to bulk as a function of cell density.

Since the cells located on the surface of the tissue have less neighbours, and thus less cooperative factors, we should expect failure to propagate from the surface towards the bulk, as observed in the experiments we report. If the radius of curvature of the shape of the tissue is much larger than the diffusion length 1/λ, a cell at the surface will be exposed to half of the cooperative factors relative to a cell in the bulk. Thus the population density near the surface n(t) can be found by simply letting β→β/2 in identities (2.8) and (2.10).

We will now estimate the velocity v=dx(t)/dt of failure propagation from the surface of a tissue of thickness (or radius) L when L≫1/λ. At t=0 the cells close to the surface have a lifetime of . Once this initial layer dies, we have a new boundary and the same process repeats itself with a new and lesser N0=N(τ). This process can be iterated in continuous time, up until time t at which the bulk collapses, or the implosion is complete. i.e. whichever is met first: t=τ or x(t)=L. We obtain the death propagation velocity as v≡dx/dt≈(2λ)−1/τ(n). That is, where b(t)=An(t)/(Φ0λ2D)∼AN0(B−Ct)1//(Φ0λ2D) is defined analogously with β.

Equation (2.12) can be integrated numerically to yield the failure depth as a function of time. Unsurprisingly, the velocity is nonzero at all times; however inspecting (2.12) reveals a non-trivial density dependence: at high cell densities the velocity is constant, independent of n. As n decreases v decreases, and reaches a minimum value. As n continues to decrease, v increases again, and diverges at some critical density that makes the parenthesis zero (figure 4). This is the critical density at which bulk death dominates surface death.

Organs have different shapes, and parts of a given organ can have varying degrees of curvature. In order to determine how the shape and curvature of a tissue affects its lifetime, we studied three idealized geometries. A tetrahedron, a cube and a sphere were formed with the same initial cell density and total cell number. The population dynamics of these three geometries were then simulated. The survival curves for these geometries can be seen in figure 5.

Figure 5.

Population dynamics of different tissue geometries. All shapes had identical volumes and initial cell densities. We used simulation parameters k=5, Φ0=175 and λ=2.

Since sharp corners are more likely to decay due to the lack of neighbours, the ‘live surface’ becomes progressively rounder as the failure propagates inwards. In the mean time, cells in the bulk also thin down and, depending on the simulation parameters, may or may not reach the critical density before boundaries collapse inward.

Discussion

Ageing is often attributed to microscopic mechanisms that cause failure at the cellular level. Our analytical, computational and experimental results support the view that ageing is due to the failure of intercellular processes as much as it is due to intracellular ones.

We have studied the details of how interactions between cells lead to the failure of the tissue. We have observed, through experiments and simulations, that damage starts from the boundaries of the tissue, and propagates inwards. At edges, external planes and points, cells have lesser neighbours compared to the bulk, and thus receive less cooperative factors. Since the number of factors received by a cell is proportional to its probability of survival, failure originates near the boundaries and propagates inwards. In addition, we found that the inward collapse is the fastest for geometries with the largest curvature or corners with the smallest angle. This trend was demonstrated experimentally and theoretically.

We have also seen that the larger death rate at the corners and edges causes the live cell surface to round up. The rounder boundary, in turn, decelerates the curvature and, thus, the local death rate. Experimentally a qualitative curvature/angle dependence of the collapse rate was demonstrated for the three corners of a triangular prism.

The most important factor determining the qualitative properties of the population collapse is the Hill parameter k. A high value of k causes an abrupt failure, while a low value of k causes a gradual loss of cells. The exact value of k should depend on the microscopic mechanism by which the cooperative factors react with the damaging agent(s). Typically, large molecules are often secreted in low numbers but have large impacts, while small molecules individually have a lesser impact, but can be synthesized quickly. If, for example, a large molecule is secreted which helps reduce radicals or peroxides, which have been shown to damage cells and their DNA [19], a single cooperative molecule may react with dozens of hazardous agents. On the other hand, if the cooperative factor is a small molecule that chemically binds to the cell, or a molecule within, which causes genetic expression or the activation of different biochemical pathways, a large number of molecules would likely be needed to increase the cell’s chance of survival. Therefore one should expect k and ϕ0 to not necessarily be independent of each other.

Our experiments involved synthetic mammalian tissues with only one cell type, and motivated the assumptions that formed the basis of our theoretical analysis. However our conclusions may also hold true for the ageing dynamics of a larger class of cooperative multicellular systems. Such systems may potentially include tumours, biofilms, colonies and microbial consortia. Thus, we reiterate that ageing is not a property of the individual components of a system, but an emergent characteristic of a strongly interacting, interdependent ensemble of components [6]. However, we should expect important qualitative differences in systems where a limiting cooperative factor is produced centrally and/or transported by a non-diffusive mechanism.

We should caution that there may be multiple cooperative factors with different diffusion lengths and different influences. Furthermore, the survival of cells may depend on non-trivial combinations of cooperative factors.

We expect multiple cooperative factors with similar diffusion lengths and influence not to change the mathematical structure of our theory, since they can be labelled as one, and summed into ϕ. If the cells vitally depend on multiple cooperative factors with different diffusion lengths and production rates, one might still get away with including in ϕ only the limiting one. However if there is a non-trivial dependence on substitutable combinations, e.g. if survival depends on (1 AND 2) or (3 AND 4), then equation (2.2) should be modified appropriately, as where the θs are Hill functions and concentrations of the four diffusive factors ϕ1,ϕ2,ϕ3,ϕ4 to be solved from equation (2.1), with different D, γ and A constants, but should otherwise have the same functional form (2.3) each. We have not studied the implications of such complications, since we have not seen any experimental evidence for them. However, our approach can be adapted accordingly, if such a complex arrangement of cooperative factors are observed to be a requirement for cell survival.

Material and methods

Isolation and maintenance of cardiac fibroblasts

Neonatal rat cardiac fibroblasts (CFs) were isolated from 2-day-old Sprague-Dawley rats (Charles River Laboratories). The rats were sacrificed by decapitation after CO2 treatment, and the hearts were immediately excised following the Institutional Animal Care and Use Committee (IACUC) guidelines of the University of Notre Dame. The excised hearts were digested in trypsin (Life Technologies) at 4°C for 16 h with gentle agitation. Next, the extracellular matrix of the hearts was further digested using collagenase type II (Worthington-Biochem) at 37°C with several cycles of agitation and subsequent trituration. Then the mixture was filtered to separate the undigested tissue pieces and the filtrate which contained cardiomyocytes and CFs was seeded into a flask. Using the differential attachment of the two cell types, CFs were separated from cardiomyocytes at the end of a 2 h incubation. These CFs were passage 1 (P1) and they were maintained by media changes every 3 days. CF cultures were passaged at approximately 80% confluency using trypsin–EDTA (0.05%; Life Technologies) and maintained in standard culture media (DMEM (Dulbecco’s modified Eagle’s media) supplemented with 10% fetal bovine serum (Hyclone) and 1% penicillin/streptomycin (Corning)).

Preparation of cell-laden hydrogels and determining cell survival

Passage 4 CFs were trypsinized and 1×105 cells were mixed 1 : 1 with 20% (w/v) PEG (Jenkem)–RGD (Bachem) conjugated PEG (PEG-RGD; PEG : PEG/RGD, 17 : 3) which contained 0.05% (w/v in phosphate buffered saline) final volume of photoinitiator (Irgacure 2959, BASF). Then, 10 μl of the mixture was sandwiched between 125 μm thick spacers, and exposed to 6.9 mW cm−2 UV irradiation for 20 s. Triangle-shaped photomasks were used during UV exposure to control the shape of the hydrogels. These conditions allow us to model three-dimensional tissues where there exist no diffusion limitations for nutrients. The hydrogels were exposed to stress by changing their media to standard culture media supplemented with 0.2 mM H2O2. Cell survival was determined using a Live/Dead assay (Life Technologies). The hydrogels were stained using ethidium homodimer-1 (stain for dead cells) every 24 h and imaged using an inverted epifluorescence microscope (Zeiss, Hamamatsu ORCA flash 4.0). For experiments intended to investigate the propagation of the decay front, the samples were stained at 0, 8, 16, 24 and 48 h with both calcein-AM (stain for live cells) and ethidium homodimer-1 and imaged using an inverted epifluorescence microscope (Zeiss, Hamamatsu ORCA flash 4.0). The number of dead cells was determined using imageJ software and the decay penetration was determined by constructing a dead cell density heat map of contrast adjusted images in Mathematica, where depth was measured as the straight-line distance from the tip to the collapse of the decay wave front. As the propagation front has a nonlinear dependence on curvature as a function of time, the straight line distance was used to create a measurement that could be compared across samples and conditions.

We performed all cell culture experiments in triplicate. Data points represent the mean value, and error bars represent the standard deviation.