Extension of Murray’s law including nonlinear mechanics of a composite artery wall

A goal function approach is used to derive an extension of Murray’s law that includes effects of nonlinear mechanics of the artery wall. The artery is modeled as a thin-walled tube composed of different species of nonlinear elastic materials that deform together. These materials grow and remodel in a process that is governed by a target state defined by a homeostatic radius and a homeostatic material composition. Following Murray’s original idea, this target state is defined by a principle of minimum work. We take this work to include that of pumping and maintaining blood, as well as maintaining the materials of the artery wall. The minimization is performed under a constraint imposed by mechanical equilibrium. We derive a condition for the existence of a cost-optimal homeostatic state. We also conduct parametric studies using this novel theoretical frame to investigate how the cost-optimal radius and composition of the artery wall depend on flow rate, blood pressure, and elastin content

. The inferred target radius of an artery then becomes a function of the flow conditions within the blood vessel (Murray 1926). In the present work, we extend this original idea of Murray so that also the material composition and wall thickness of the artery are determined by a minimum work principle. This development ties together with the previous work (Satha et al. 2014), where we studied how local changes in volumetric blood flow or pressure, due to, for instance, disease, injury, and surgery, trigger growth and remodeling (Humphrey 2002) toward a homeostatic target state. In this paper, we develop a theory that determines this target homeostatic radius, wall thickness, and material composition, the artery wall being a composite of different constituents with nonlinear material properties (Holzapfel et al. 2000). In order to keep the theory as simple as possible, we assume the vessel to be of cylindrical shape, and we use a theory for thin-walled structures.
The blood vessel wall mainly consists of elastin, collagen, and smooth muscle (Boron and Boulpaep 2008, pp. 473-481). Thus, we model the vessel wall as a composite of multiple orthotropic, nonlinear elastic materials that deform together as the vessel stretches in the circumferential direction due to the transmural pressure, as described in the literature (Humphrey and Rajagopal 2002;Gleason and Humphrey 2004;Satha et al. 2014). The target composition and radius are assumed to minimize the cost-that is, the power per unit length of blood vessel-required to maintain and pump the blood contained within the vessel and to maintain the materials of the vessel wall, as previously proposed (Taber 1998;Klarbring et al. 2003;Liu and Kassab 2007). The goal function of the system is then taken to be this cost function subject to the constraints imposed by the mechanical equilibrium of the vessel wall. Since the elastin content changes very slowly in the vascular system of adult individuals (Tsamis et al. 2013), the amount of elastin is essentially beyond the control of the growth and remodeling process. Therefore, we regard the amount of elastin as a parameter to the system. The goal function is then parameterized by the blood pressure, the volumetric flow rate, and the amount of elastin. These parameters, in turn, are functions of time, and their fluctuations lead to fluctuations of the target geometry and composition.
Experimental studies show that an increased blood pressure p increases the thickness of the vessel wall through growth and that the vessel adapts to achieve a homeostatic state (Matsumoto and Hayashi 1996;Hu et al. 2007). These studies also show that changes in blood pressure affect the material composition of the vessel wall. Similarly, the volumetric flow rate u has a strong impact on a blood vessel's radius and composition: The radius r is increased when the flow rate is increased, so that the shear stress of the fluid on the epithelial cells, that is, the interior lining of the vessel wall, is kept at a homeostatic state (Brownlee and Langille 1991). On a longer timescale, the material composition of the vessel wall also changes with increased flow rate (Kubis et al. 2001). It was suggested in an early work by Murray (1926) that the target dimensions of the blood vessel are governed by the minimization of metabolic power needed to maintain the materials of the vascular system and to overcome the hydrodynamic resistance from the vessel for a given demand of supplied blood. This minimization principle leads to Murray's law which is in fair agreement with the experimental data (Sherman 1981;Taber et al. 2001). Later, Murray's law was modified by taking the metabolic cost of the vessel wall into account (Taber 1998), including the active behavior of smooth muscle. This latter approach relates the shear stress of the homeostatic state to the pressure, the thickness of the vessel wall, and the degree of smooth muscle metabolism. Klarbring et al. (2003) and Liu and Kassab (2007) have further developed the cost function approach by considering minimization of the cost for the vascular tree as a whole in their formulations.
To the knowledge of the authors, the fact that the artery wall is composed of several constituents with orthotropic, nonlinear properties (Holzapfel et al. 2000) has not been considered in previous studies of the cost-optimal geometry and composition of artery walls. Because the elastin content of the artery is essentially unchanging at the timescales of growth and remodeling (Tsamis et al. 2013), there is not a unique optimal target composition of the artery wall for a given set of flow parameters; the optimal state depends on the given amount of elastin, and its slow variations due to degradation. The target composition may then be coupled to the material properties of the composite artery wall.
To find the cost-optimal geometry and composition of an artery with a nonlinear mechanical behavior, it is necessary to consider a mechanical model of the artery wall in conjunction with a cost function derived from the power required to maintain the materials and blood flow of the artery. We briefly outline the mechanical model, based on constrained mixture theory (Humphrey and Rajagopal 2002;Gleason and Humphrey 2004;Satha et al. 2014), in Sect. 2.1. This yields an equilibrium equation that relates the transmural pressure to the vessel geometry and composition of a homeostatic state. A description of the principle of cost-optimization for the artery wall follows in Sect. 2.3, and a goal function is subsequently formulated, whose minima correspond to a minimal cost of homeostatic states that satisfy the equilibrium equation (Sect. 2.4). We analyze how the cost-optimal state of the vessel varies with volumetric flow rate, pressure, and elastin content in Sect. 3.

Constrained mixture thin-walled tube theory
We consider a cylindrical tube composed of a mixture of n materials, whose respective mechanical properties are represented by their strain energy functions ψ k , k = 1 . . . n. A constrained mixture theory is used, implying that all constituents have the same deformation. This deformation, with respect to a given, fixed reference configuration, is represented by a circumferential strain λ and a supposed constant axial strain δ. For a pressure difference p between the interior and exterior of the tube, integration of the standard radial equilibrium equation gives where ρ is the radial coordinate which varies between an inner radius ρ 0 and an outer radius ρ 1 . For an incompressible material, the stress difference between circumferential stress σ ϕ and the radial stress σ ρ can, cf. Holzapfel and Ogden (2003), be written where φ k denotes the volume fraction of constituent k. Introducing Eq. (3) into Eq.
(2) and making a thin-walled tube assumption, cf. Satha et al. (2014) for details, result in where R is the radius of the, now thin-walled, reference configuration, and A k is the effective reference area obtained by multiplying the volume fraction φ k by the total reference cross-sectional area. The radius of a deformed, thin-walled tube is expressed as r = λR. Essentially following Baek et al. (2006), we take the effective areas to be represented by where A k (0) is the original effective area of constituent k, Q k (t) is the fraction of constituent k that was produced before time 0 and remains at time t, A k (t) ≥ 0 is the rate of production of effective area at time t, and q k (t) ≥ 0 is a monotonically decreasing survival function such that q(0) = 1. By assuming that materials created at different time instances contribute to the strain energy in proportion to the remaining area fractions, we obtain (Baek et al. 2006) where, Ψ k λ k (t, τ ) is the strain energy density with respect to a natural, stress-free configuration and characterizes the nonlinear, elastic behavior (Baek et al. 2006). Also, λ k (t, τ ) is the stretch at time t for materials produced at time τ . Hence, (Baek et al. 2006) The ratio λ(t)/λ(τ ) is the stretch developed during the time interval [τ, t], and G k h is the homeostatic prestretch of constituent k, which means the material may attain prestretch at the time of production.

Timescales and homeostatic conditions
We recognize different timescales in the process of growth and remodeling of the vascular system. The high-frequency scale is that of the heartbeat. It was shown in Satha et al. (2014) that Eq. (4) is approximately valid for average quantities if the change of A k is taken to be much slower than that of the heartbeat timescale. Moreover, we distinguish between two processes in the slow change in A k . First, there is the change of homeostatic values. Secondly, there is the process of approaching these homeostatic target values when, say, a perturbation of the state occurs. The stability of the second type of process was previously investigated in Satha et al. (2014). Complementary to this, in the present paper, we study the target homeostatic state and its dependence on the imposed flow conditions. Such states are defined by a time-constant stretch λ(τ ) ≡λ as well as a time-constant composition of materials A k ≡Â k . There are two classes of constituents for which steady-state conditions are possible (Satha et al. 2014): The set of constituent indices belonging to class (i) and (ii) are denoted by S i and S ii , respectively. Equations (5) and (6) result in (Satha et al. 2014) Here, G k Introducing Eq. (8) into a time-averaged version of Eq. (4), and evaluating for λ = λ(t) = λ(τ ) =λ and for A k = A k (t) =Â k , we get (Satha et al. 2014) where is called the homeostatic stress. Note that the homeostatic state is associated with a constant homeostatic stress for materials with a finite turnover. Here and in the following, we use the notation d f (s) = d f/ds and d n f (s) = d n f /ds n .

Principle of cost-optimization
As proposed by Murray (1926), it is assumed herein that the blood vessel growth and remodeling strive toward costoptimization of the vascular system. This assumption has been widely used in previous modeling work (Taber 1998;Klarbring et al. 2003;Liu and Kassab 2007). In this work, we take the target homeostatic state to be governed by such an optimization rule. We assume that the metabolic cost of the materials that constitute the vessel wall is proportional to the amount of each constituent, i.e., there are constants α k such that this cost per unit length of blood vessel in the homeostatic state can be written with the units of power per unit length. Since the homeostatic stress of smooth muscle is constant (Sect. 2.2), it is possible to represent the stress-dependent upkeep of smooth muscle (Taber 1998;Liu and Kassab 2007) by the constant α k . There is also a metabolic cost for the blood. This is again taken as proportional to the volume, i.e., it is proportional to πr 2 δ. Since r =λR, and since a constant axial stretch δ is considered, there is a constant β such that the metabolic cost of the blood per unit length of the blood vessel can be written We have β = πδα b where α b is the metabolic power per unit volume of blood. Finally, we take into account the energy per unit time consumed by the heart to maintain a certain volumetric flow rate. If we assume that the Hagen-Poiseuille equation governs the flow, the power per unit length of blood vessel required to overcome the viscous drag is (Taber 1998) where u is the volumetric flow rate, and η is the dynamic viscosity of the blood, which is assumed to be a Newtonian fluid. There is thus a constant γ = 8η/π such that the cost is per unit length of blood vessel. The total cost P per unit time and length is obtained as the sum of these contributions, becoming

The optimization problem and its minima
The problem we are considering is thus to minimize the total cost P under the constraint that the equilibrium condition, Eq. (9), is satisfied. This problem can be rewritten as an unconstrained optimization problem by taking an arbitrary j ∈ S i and rewriting Eq. (9) asÂ Thus,Â j =Â j (λ,Â 1 , . . . ,Â j−1 ,Â j+1 , . . . ,Â n ), and when substituted into the expression for P in Eq. (15), we get the goal function The target homeostatic state is now given by the unconstrained minimum of f , assuming that this minimum occurs for positive values of all variables. The model is next simplified by assuming that the blood vessel wall consists of two constituents only: elastin, k = 'e' ∈ S ii , and components with a finite turnover including collagen and smooth muscle, k = 't' ∈ S i . This classification incorporates the assumption that the elastin content is essentially constant over time (Tsamis et al. 2013), while other constituents have a substantially faster turnover, with a timescale of approximately 2 months (Nissen et al. 1978;Martufi and Gasser 2012). Smooth muscle is metabolically more expensive than collagen, and it is present in the vascular system to help pumping blood and to control high-frequency adaptation to changing demands of blood. The fraction of smooth muscle is then likely related to the fluctuations of the flow conditions rather than their time-averaged values. However, these dynamics are beyond the scope of this study, and we introduce the simplifying assumption that the ratio of the amount of collagen to the amount of smooth muscle is constant for any given artery.
For the two constituents, 'e' and 't,' we can express the equilibrium equation (16) aŝ where σ t h is a constant homeostatic stress, and Eq. (10) was used to express σ e h (λ). Substituting Eq. (18) into the total cost P gives the goal function This cost function, retaining only nonconstant terms, becomes and the gradient of the goal function is Straight-forward differentiation of Eqs. (20) and (18) yields The optimal target homeostatic composition of a blood vessel is found at a stationary minimum point defined by Using that ∂ P/∂Â t = α t is constant, the second derivative of f is We note that Then, d 2 f /dλ 2 > 0 when α t = 0. If α t > 0, we must consider the sign and magnitude of d 2Ât /dλ 2 : Whether or not this expression is positive at a stationary point can be evaluated when the material model is instantiated. This will be done in Sect. 3.1. However, qualitative insight can be gained by equivalently writing Eq. (28) as Thus, in case the elastin stress σ e h is proportional toλ, so that the second term vanishes, the stationary point will always be a minimum point. On the other hand, if the elastin has a strainstiffening behavior, then d 2Ât /dλ 2 may become negative. Particularly, this may be the case for small pressures.
If we assume that the metabolic cost of the vessel wall is much smaller than that of the blood, α t ≈ 0. Then, d f/dλ = 0 gives consistent with Murray's law (Murray 1926). This result can be inserted into Eq. (18) to give a closed expression for the optimum amount of materials with finite turnover. For a finite metabolic cost of the vessel wall, α t > 0, the stretch at the stationary point of the goal function must be computed numerically for any nontrivial choice of strain energy function Ψ e .

Results and discussion
The cost-optimal target geometry and composition of the vessel wall are found at the minimum stationary point of the goal function. The locus of this stationary point depends on the parameters of the goal function, including pressure p, volumetric flow rate u, elastin contentÂ e , and parameters related to the material model for elastin. These parameters vary within a population as well as with time for each individual due to, e.g., aging, changes in body mass, medical treatments, or the development of diseases. In Sect. 3.2, we perform parameter studies to quantify these variations in the optimal state. However, we first need to be explicit about the material model and its parameters.

Parameter identification and material model
The parameters of our model are quantified using data for the radial artery (arteria radialis) and the common carotid artery (arteria carotis communis). Previous in vivo measurements on normotensive subjects are used, giving ensemble averages for the vessel radiusr , total areaĀ of the cross section, average blood pressurep, and volumetric flow ratē u, as compiled in Table 1. The composition, described by the fraction of elastin φ e and the fraction of other materials φ t , is estimated using histological data from the literature, as described by Satha et al. (2014). We use histological data from Li et al. (2008) for the radial artery and from Sommer et al. (2010) for the carotid artery (Table 1). The stretches in the circumferential, radial and longitudinal directions are λ k , (λ k δ) −1 and δ, yielding the Cauchy-Green tensor (Holzapfel et al. 2000; Ogden 2010) with invariants I k 0 = tr C k and I k 1 = (λ k ) 2 . As previously described (Satha et al. 2014), the strain energy density of the elastin fraction is taken to be isotropic (Holzapfel and Ogden 2010):  Giannattasio et al. (2001) b Laurent et al. (1994) c Estimated by Satha et al. (2014) using histology from Li et al. (2008) d Numerical fit by Satha et al. (2014) to data from Laurent et al. (1994) and Girerd et al. (1998) e  f Likittanasombut et al. (2006) g Bussy et al. (2000) h Bussy et al. (2000) with a correction for a misrepresented unit i Sommer et al. (2010) j Numerical fit using the method of Satha et al. (2014) with data from Bussy et al. (2000) while the strain energy density of the composite of other constituents is taken to be orthotropic (Holzapfel and Ogden 2010): where c 1 > 0 Pa is a constant and c 2 > 0 is a nondimensional constant. Parameter identification for the radial artery was performed in a previous study (Satha et al. 2014) by leastsquares fitting the two-constituent material model to experimental data (Laurent et al. 1994;Girerd et al. 1998), giving the material parameters shown in Table 1. Using Eq. (10), these parameters yield σ t h = G t h dΨ t (G t h ) = 38.1 kPa. The fitting procedure described by (Satha et al. 2014) is used herein to obtain the parameters of the carotid artery from the data of Bussy et al. (2000), with the Young's modulus of the unloaded wall of the carotid artery estimated to 0.3 MPa, similar to the value for the brachial artery (Kinlay et al. 2001). The resulting material parameters for the carotid artery are compiled in Table 1 and give σ t h = G t h dΨ t (G t h ) = 46.3 kPa. We also choose the constant longitudinal stretch to be δ = 1.
The parameters, α t , β, and γ , of the goal function are obtained from the literature. Liu et al. (2012) estimate α b = 51.7 W/m 3 for human blood, giving β = 0.16 kW/m 3 . With a Newtonian fluid assumption, the dynamic viscosity of human blood at 40 % hematocrit is η = 3.2 mPa·s (Boron and Boulpaep 2008), giving γ = 8.1·10 −3 Js/m 3 . The metabolic coefficient α t is assumed to be dominated by smooth muscle and has an active and a passive component, with the active component proportional to the stress of that constituent (Taber 1998). We thus write where α w and k w denote the passive and active metabolic constants, respectively. These constants were estimated by Taber (1998) to be α w = 764 W/m 3 and k w = 0.00872 s −1 for the porcine carotid artery, giving α t = 1.1 kW/m 3 and α t = 1.2 kW/m 3 for the radial and carotid artery, respectively. We take these values for α t as order of magnitude estimates for human arteries and investigate different values α t = {0.0, 0.1, 1.0} kW/m 3 in the parametric studies below.

Parametric studies
In this section, we consider the effects of the volumetric flow rate, pressure, and elastin content on the radius r of the blood vessel and on the amount of constituentsÂ t with a finite turnover. The parameter α t , controlling the cost of the 't'type wall materials, is varied to highlight its effect on the vessel dimensions and composition. The target state for each set of parameters is found numerically by solving d f/dλ = 0 forλ using Eqs. (19) through (24), and then computingÂ t using Eq. (18). From the point of view of growth stability, it is of great interest to assess whether the stationary points of the goal function are minima. With the prototypical values from Table 1, we have evaluated Eq. (28) for a wide range of the radius 0.1r < r < 3r and the pressure 0.2p < p < 3p and found that d 2Ât /dλ 2 > 0 within these ranges for both the radial and the carotid arteries. This means that the second derivative of the goal function with respect toÂ t is strictly positive, asserting that the corresponding stationary points are indeed minima. Note that this validation was conducted for one particular choice of material model. An enhanced strain-stiffening, e.g., due to an anisotropic elastin fraction, would lead to greater nonlinearity in the strain energy density which would threaten the existence of the minimum. Therefore, we cannot exclude that there exists some physiological conditions at which the minimum of the goal function is lost.

Volumetric flow rate
It has been established experimentally that the volumetric flow rate has a strong impact on the blood vessel radius (Brownlee and Langille 1991;Kubis et al. 2001) and com-  (Kubis et al. 2001). In our theoretical framework, this is manifested as a flow rate dependence of the stationary point of the goal function. The vessel radius r and the amount of composite materialsÂ t are plotted as functions of u in Fig. 1a, b (radial artery) and Fig. 1c, d (carotid artery) for different values of α t = {0.0, 0.1, 1.0} kW/m 3 and a constant pressure p =p.
When the cost of wall materials is taken to be zero, α t = 0, the variations of r with u follow Murray's law, r ∝ u 1/3 ( Fig. 1a-c, solid line). Murray's law overpredicts the average vessel radiusr given in Table 1 for the average flow rateū. When the wall material is assigned a finite metabolic cost, Murray's law is modified to suppress the use of wall materials and thus reduce the radius to a more realistic value (Fig. 1a-c,  dotted lines). Interestingly, this also introduces a lower bound on the vessel radius, which does not fully contract even at a vanishing flow rate.
When examining the relationÂ t (u) for the radial artery (Fig. 1b) and the carotid artery (Fig. 1d), it becomes clear that A t > 0 for all flow rates investigated. There is a minimum of A t (u) that corresponds to a zero of dÂ t /dλ. For a constant pressure, the amount of materials in the vessel wall is a rather weak function of the flow rate.
The rise in the amount of 't'-material for low-volume flows corresponds to the elastin being in a state of compression, requiring additional 't'-material, with constant stress σ t h , to balance the pressure p. However, this may be an artifact of

Pressure
When the cost of the wall materials is taken to be zero, α t = 0, and Murray's law governs the target state, the pressure does not have any effect on the vessel radius, as shown for both the radial and the carotid arteries by the solid lines in Fig. 2a-c. Also, it is observed in Fig. 2b-d thatÂ t is linear in pressure, which is trivially explained by the need to balance the pressure at a constant circumferential stress σ t h in the 't'-fraction of materials. Examining the solid lines in Fig. 2b-d, we note thatÂ t becomes negative when the pressure is sufficiently reduced. Below this limiting pressure, no realizable homeostatic state can be found which reproduces the prediction of Murray's law. This constitutes a lower limit of pressure for Murray's law. This may also be an artifact of the simplistic model for the strain energy density of elastin in compression, as discussed in Sect. 3.2.1.
Under normal circumstances, with a typical pressure p = p, assigning a finite cost to the wall material, α t > 0, leads to a more narrow blood vessel (Fig. 2a- (Table 1). A narrow blood vessel reduces the force per unit length of the vessel wall and thus allows for a thinner wall, which saves expensive materials. It is interesting that the vessel radius increases when the blood pressure is reduced: A reduced blood pressure at a sustained volumetric flow rate then reduces the mechanical stability of the vessel and increases the risk of vessel collapse. The dramatic increase in the radius at low pressure is not physiological, since it occurs at states withÂ t < 0 ( Fig. 2b-d, dashed and dotted lines), which can never be achieved.

Elastin content
Although the elastin content is essentially constant (Tsamis et al. 2013), it may degrade over very long timescales, e.g., the lifetime of an individual. This motivates a study on how variations-particularly reductions-in elastin content affect the homeostatic target vessel geometry and composition. Figure 3a, c show how the radii of the radial and carotid arteries, respectively, vary with the elastin content. When α t = 0, the vessel radius is maintained at a constant level, owing to the fact that the elastin content does not enter into Murray's law (Fig. 3a-c, solid line). Degradation of elastin is compensated for by an increase in the amount of other mate-rialsÂ t . It is shown in Fig. 3b-d thatÂ t (solid line) increases linearly whenÂ e is reduced. That is, degraded elastin is simply replaced by other materials to balance the transmural pressure.
For the case α t = 0, elastin is replaced by metabolically more expensive materials. This is predicted to lead to a reduction of the vessel radius when elastin degrades (Fig. 3a-c, dashed and dotted lines).

Comparison between radial and carotid artery
To demonstrate the general applicability of the proposed model, two types of arteries, the radial artery and the common carotid artery, are compared. These arteries are very different in terms of diameter and blood flow, but have a similar transmural pressure. The fraction of elastin is much greater in the carotid artery ( Table 1).
The predicted variation of the vessel radius r with u deviates significantly from Murray's law for the radial artery (Fig. 1a), whereas the Murray's law appears to hold much better for the carotid artery (Fig. 1c). The same conclusions can be drawn for the amount of 't'-materialsÂ t (Fig. 1b-d).
In the cases of pressure dependence and elastin content dependence, the radial and carotid arteries display the same qualitative behavior, which clearly differs from Murray's law .

Conclusions
The design of the vascular system is assumed to be governed by the physiological principle of minimum work (Murray 1926). It is thus an optimization process that governs the architecture of arteries. On this basis, we have formulated a theoretical frame that extends Murray's law to include growth and remodeling, and the nonlinear mechanics of the artery wall. A goal function, novel to this application, is formulated using an expression for the power required to pump blood and the total metabolic power needed to maintain the blood and the wall of the artery.
We have shown that there exists a minimum stationary point for a wide range of the volumetric flow rate and the pressure around the prototypical parameter values for the radial and the common carotid artery. In theory, however, this minimum could be lost for a strongly strain-stiffening elastin fraction.
Taking the cost of wall materials into account reduces the radius of the target homeostatic state and also renders this target radius pressure-dependent. A reduction in the amount of elastin in the artery wall reduces the radius of the target homeostatic state.
The greatest value of the present work may be its ability to depict the variations of the target homeostatic state under dynamic flow conditions. This theoretical frame can then be integrated into models for growth and remodeling (Satha et al. 2014;Taber 1998) to capture the coupled dynamics of remodeling and fluctuation of the target state.