S. Negin Mortazavi Department of Mechanical Engineering, University of Texas at Dallas, Richardson, TX 75080 e-mail: [email protected]

Foteini Hassiotou Assistant Professor School of Chemistry and Biochemistry, University of Western Australia, Crawley, Western Australia 6009, Australia e-mail: [email protected]

Donna Geddes Associate Professor School of Chemistry and Biochemistry, University of Western Australia, Crawley, Western Australia 6009, Australia e-mail: [email protected]

Fatemeh Hassanipour1 Assistant Professor Department of Mechanical Engineering, University of Texas at Dallas, Richardson, TX 75080 e-mail: [email protected]

1

Mathematical Modeling of Mammary Ducts in Lactating Human Females This work studies a model for milk transport through lactating human breast ducts and describes mathematically the mass transfer from alveolar sacs through the mammary ducts to the nipple. In this model, both the phenomena of diffusion in the sacs and conventional flow in ducts have been considered. The ensuing analysis reveals that there is an optimal range of bifurcation numbers leading to the easiest milk flow based on the minimum flow resistance. This model formulates certain difficult-to-measure values like diameter of the alveolar sacs and the total length of the milk path as a function of easyto-measure properties such as milk fluid properties and macroscopic measurements of the breast. Alveolar dimensions from breast tissues of six lactating women are measured and reported in this paper. The theoretically calculated alveoli diameters for optimum milk flow (as a function of bifurcation numbers) show excellent match with our biological data on alveolar dimensions. Also, the mathematical model indicates that for minimum milk flow resistance the glandular tissue must be within a short distance from the base of the nipple, an observation that matches well with the latest anatomical and physiological research. [DOI: 10.1115/1.4028967]

Introduction

Evolution of many biological flow systems has given rise to structures that provide better access paths for the flow. Specifically, tree-shaped mass transfer structures are common in biological organs, e.g., the lung, the kidney, the breast, and other vascular tissue. A deeper understanding of this complex branching architecture of organs is essential for elucidating both the physical basis of the transport processes and the pathology of many human diseases, such as agenesis and asthma in the kidney and lung, as well as breast carcinoma, a disease of terminal ductal lobular units [1]. Many aspects of the diagnosis, healing, and investigation of the origin of breast diseases would significantly benefit from a detailed knowledge of the breast branching morphogenesis [2,3]. A great number of conventional breast conditions, such as ductal blockage, breast engorgement, breast abscess, and galactocele, can render breastfeeding difficult and sometimes impossible. There are many factors that contribute to diseases of the breast, but among the most important are those that relate to the mechanical properties of the breast, e.g., ductal physical properties. For milk production in the breast (diffusion in sacs), the factors that move the fluid of milk in the ducts toward the nipple (sucking pressure) are all essential mechanical processes. Despite the clear advantage of a fundamental understanding of the branching structure properties in the human breast ductal system, this has not yet been the subject of careful mathematical modeling and evaluation. Thus, the present study produces a mathematical model that reveals the interactions between many parameters in this complex biological system, providing new and interesting insights. Many other biological organs have been subjected to extensive theoretical studies of the mechanical properties of their flow structure, e.g., lung [4–7], nasal airway [8–11], arteries in vessels [12–19], brain [20,21], and liver [22–24]. However, to the 1 Corresponding author. Manuscript received April 16, 2014; final manuscript received October 17, 2014; published online June 3, 2015. Assoc. Editor: Alison Marsden.

Journal of Biomechanical Engineering

knowledge of the authors, the mechanisms of milk flow in the ductal system of the lactating breast have not been the subject of a fundamental mathematical study. The studies of mammary duct structures have been done either through imaging (ultrasound) or visualizing (anatomization) the three-dimensional (3D) ductal structure, mostly by biologists or radiologists. The first model is produced by injection of each ductal/lobular system of women who died during lactation with colored wax from the nipple openings before dissection [25]. The lobes are identified using different colors, and ducts are isolated from each other and clearly distinguished. In another study, a computerized 3D anatomy of the ductal system has been developed by manually tracing the ducts from histological sections onto acetate sheets and by analyzing and reconstructing the breast ducts [26]. Also a similar computer-assisted 3D construction of the ductal system of a postmenopausal breast is developed for an entire breast, resected by subcutaneous mastectomy from a 69-year-old woman, with ductal carcinoma [27]. The number, distribution, and anatomic properties of the ductal systems of the breast (from the nipple orifices to the terminal lobular units) have been determined using six complementary in vivo and in vitro approaches [3]. Using ultrasound imaging, the anatomy of the lactating breast with the aim of elucidating clinical implications particularly for breast surgery has been investigated [28]. A 3D phantom model with realistic and randomized anatomical features is developed to generate a realistic model for a breast imaging system [29]. The mathematical model proposed in this paper can be utilized to deduce the milk flow from the mechanical properties of the bifurcated ducting structure in the lactating breast. Our model represents the total path from the production zone (alveoli) to the outlet zone (nipple), leading to closed-form expressions that yield the milk flow parameters from the geometric parameters of the model. An immediate example of the applications of this model is the calculation of flow resistance as the ratio of the pressure drop and the milk flow rate. Flow resistance is an important mechanical property of the breast ductal system that must be successfully overcome by the sucking pressure. Calculation of milk-way resistance

C 2015 by ASME Copyright V

JULY 2015, Vol. 137 / 071009-1

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

is made possible by the notion that the ductal tree behaves like a connected conduit through which flows an incompressible fluid (milk). This paper finds a range for the optimum number of branchings based on the easiest milk flow through ducts, in other words, the minimum overall flow resistance. It is demonstrated that the combination of convective and diffusive mechanism allows for the maximization of milk transfer through the mammary ducts in order to achieve the best use of this system.

2

Breast Development

Human breast is a dynamic organ that does not reach its full maturity unless a woman experiences pregnancy and childbirth [30]. During pregnancy, the mammary glands and the terminal ductal lobular units significantly expand in anticipation of lactation. The completed milk flow system in the lactating breast looks like a tree. The leaves of the tree are the alveoli that form grapelike spherical cavity structure (Fig. 1). Alveoli generate milk from the nutrients and water that are diffused from the adjacent blood stream. A cluster of alveoli is called a lobule. Approximately, 20–40 lobules combine to make a lobe [31]. Several studies [26,27,32] indicate that more than 90% of nipples have 15–20 lobes, out of which only 5–9 are true mammary duct orifices and the rest are sebaceous glands which vary in length from 1 cm to 4 cm with no branches and connections to the ducts [3,28]. The peripheral ducts cover the central ducts in a radial fashion in a way that these two systems do not connect to each other. Milk is distributed to different parts of the nipple from the alveoli through the ducts. Total ductal cross-sectional surface area decreases from alveoli to nipple; therefore, the milk flow accelerates to reach the nipple. There are millions of alveoli in the lactating breast. This complicated structure with various duct lengths, duct diameters, and other physical properties is required for the delivery of milk to the sucking infant in a trouble-free and efficient manner. The present work studies the lactating breast medium from a macroscale point of view as a continuum which includes solid phase (tissue) and fluid phase (milk) in alveoli and mammary ducts.

3

Theoretical Modeling of Flow Resistance

Milk flow in the lactating breast has two paths: diffusion of blood-derived components (protein, glucose, and fat) through alveoli (diffusive resistance) and flow through mammary ducts (convective resistance). The total minimum flow resistance is obtained when the structure is arranged in a way that the flows with high resistance (alveoli) occupy the smallest scales of the flow system, while the flows with the lower resistance (ducts) inhabit the larger scales. To have a better insight about milk flow in the lactating breast, each of the two resistances are discussed and evaluated in Secs. 3.1 and 3.2. 3.1 Convective Resistance. Milk flows along a conduit with circular cross section (Fig. 2). Considering the largest values of initial duct diameter (2 mm) [28], maximum observed milk flow rate (4.8 gr/min) [33–35], and the minimum measured milk viscosity (1.66  106 m2/s) [36], the maximum Reynolds number of the milk flow in the conduits will be around 30 which can be considered as laminar flow. Therefore, the milk flow resistance of each of the two ducts arising from bifurcation (n) can be represented by the well-known Hagen–Poiseuille equation [37] Rd ðnÞ ¼

Anatomy of the lactating breast

071009-2 / Vol. 137, JULY 2015

(1)

where (Rd) is the resistance of a single duct, (DPn), ðm_ n Þ; ðqÞ; ðlÞ; ðDn Þ, and (Ln) stand for pressure drop, mass flow rate, milk density, dynamic viscosity, duct diameter, and duct length (in generation (n)), respectively. This formula often appears in texts of pulmonary physiology [38]. The convective resistance calculation is based on a single suckling pulse. In a separate work, the authors have studied the behavior of milk flow in mammary ducts with pulsating flow during the total period of breastfeeding [35]. For the geometry of the model, a simple structure of the mammary ductal tree has been considered with dichotomous branching scheme for each lobule (Fig. 3(left)). Each parent duct is divided into two identical daughters and this division is repeated many times, with each successive division being called a generation. Assuming complete symmetry (Fig. 3(left)), the value of two dimensionless ratios, (Dn/Dn–1) and (Ln/Ln–1), which are the ratio between consecutive duct diameters and duct lengths, can be found by constructal law which has been used in several other natural and biological domains to simplify the theoretical calculations of a complicated branching system [39–43] Dn ¼ 21=3 Dn1

Fig. 1

DPn 128lLn ¼ m_ n pqD4n

and

Ln ¼ 21=3 Ln1

(2)

Fig. 2 Cross-sectional view of human lactating breast duct

Transactions of the ASME

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

Considering the fact that the main flow stream ðm_ n Þ is divided into two separate streams ðm_ nþ1 Þ, where ðm_ nþ1 ¼ m_ n =2Þ, and combining Eqs. (1), (2), and (6), the variation of kinetic energy per unit mass in generation (n) would be Den ¼

pD40 22n=3 ðDPÞ2 218 l2 L20 q

(7)

Rearranging Eq. (7) derives the bifurcation resistance in the joints in generation (k) as Rbf ðkÞ ¼

Fig. 3 Dichotomous branching of ducts and representative model

where (Dn) and (Ln) are, respectively, ductal diameter and length in generation (n). With the above assumption, the number of milk paths in generation (n) is 2n and the total convective resistance can be obtained by adding the individual resistances in series or in parallel (Fig. 3(right)). Assuming the same pressure drop across each of the two daughters in the same generation, the convective resistance to the laminar flow in each single duct in generation (0) to (n) is Rd ð0Þ ¼ 20

Rgen bf ðkÞ ¼

(3)

Rtot bf ¼

k¼0

n1 X

Rgen bf ðkÞ

¼

n1 X

m_ 0 22k=3 2 D4 16pq 0 k¼0

m_ 0 ð1  22n=3 Þ ! 1010 755qL0 l

(10)

Therefore, the total convection resistance (ductal þ bifurcation) for the mammary gland with nipple at (n ¼ 0) and considering (n  1) bifurcation levels is given by

The generations act in series, therefore, the overall resistance increases with the number of generations as follows: ¼

(9)

Similar to the series–parallel approach for the total ductal resistance, the total bifurcation resistance will be

¼

where (D0) and (L0) are the diameter and the length of the first mammary duct, respectively. In each generation (k), we note that due to symmetry the pressure at the upstream ends of all the ducts in the same generation is equal, and similarly, the pressure at the downstream end of these ducts is also equivalent. Therefore, for computational purposes, one may consider that these ducts operate in parallel. Thus, using the physical laws of combining parallel resistors, the total resistance contributed by generation (k) is calculated Rk (4) Rgen d ðkÞ ¼ k ¼ Rd ð0Þ 2

Rgen d ðkÞ

m_ 0 Rbf ðkÞ ¼ 22k=3 2k 16pq2 D40

k¼0

128lL0 pqD40

Rd ð2Þ ¼ 22 Rd ð0Þ .. . Rd ðnÞ ¼ 2n Rd ð0Þ

n X

(8)

Thus, the total bifurcation resistance in generation (k) of bifurcation, by taking into account that there are (2k) mammary ducts in generation (k), can be obtained as

Rd ð1Þ ¼ 21 Rd ð0Þ

Rtot d

m_ 0 DPk Dek ¼ 2k=3 ¼q m_ n k m_ k 16pq2 D40

128lL0 ¼ ðn þ 1ÞRd ð0Þ ¼ ðn þ 1Þ pqD40

Rtot cn ¼ ðn þ 1ÞRd ð0Þ

3.2 Mass Diffusion Resistance. During lactation, the infant’s sucking stimulates the hormone Oxytocin and prompted by that, alveoli walls start contracting [44–46] and the bloodstreams (milk ingredients) approach the capillaries near the breast tissue into alveoli to form the building blocks of milk. The bloodstream moves through the cells that line the alveoli and into the sacs [47,48]. Assuming a steady diffusion of fluid with density, (q), through the alveolous’ wall with thickness of (d), and surface area _ (temporal rate of arrival of of (A), the rate of the diffusion ðmÞ fluid mass in the sac) in one dimensional form will be [49] m_ ¼

(5)

(11)

qDf A DM d

(12)

where (u) is the flow velocity profile in ducts.

where (Df) is the diffusion coefficient (m2/s) [50] of blood through the membrane walls. (DM) is deriving potential for flow of mass through the alveoli’s wall [49]. Considering the fact that the alveoli muscle contraction, (Dpsac), leads to flow of bloodstream and as a result of alveoli volume changes (expansion), the deriving potential (DM) can be calculated as (Dpsac/E) [49] where (E) is the bulk modulus of elasticity of the alveoli tissue. (d) represents the total distance from blood side to the milk side in alveoli (Fig. 4). (d) includes the capillary cell lining, the basement membrane, and the two epithelial layers. Recent studies [51] have shown that the thickness of alveoli is approximately d ¼ 0.2dsac. Alveolar sacs can have various shapes, but in a simple approach they can be considered approximately spherical with the alveolar diameter set as (dsac). Therefore, the total mass of blood components diffusing to the alveolar sacs will be

Journal of Biomechanical Engineering

JULY 2015, Vol. 137 / 071009-3

Bifurcated junctions also add extra resistance to the convective milk flow system; however, this value is in the order of 1010 and can be neglected compared to the ductal resistance (Eq. (5)). This specific resistance can be calculated by determining the energy dissipation and the pressure drop occurring between the upstream and the downstream of the bifurcated joint [42]. Applying the variation of kinetic energy per unit mass (e) for two inlet and one outlet in a joint, we will have [42] Den ¼

q 2

ð

u3 dA  _ m nþ1 nþ1

ð

u3 dA m n _n

 (6)

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

Fig. 4 The architecture of each lobe, consisting of several lobules

m_ ¼

pdsac qDpsac Df 0:2E

(13)

Our model consists of (2n) sacs in each lobule, and (a) lobules per lobe (Fig. 4); therefore, the total number of sacs per lobe is (2n  a) resulting in a proportional mass flow rate (Eq. (14)). According to the latest literature on human breast anatomy [31], the number of lobules per lobe is estimated to be between 20 and 40. Therefore, the total diffusive mass flow rate as a function of pressure drop will be n X

m_ sac ¼ ð2n  aÞ

0

pdsac qDpsac Df 0:2E

(14)

and the diffusive resistance will be Rdf ¼

0:2E ð2n  aÞpdsac qDf

Therefore, the diffusion resistance for one lobule in Eq. (15) will be   0:0412E ð2nþ1Þ 2 3 (18) Rdf ¼ apL0 qDf

4

Results and Discussion

The convective resistance is an increasing linear function of (n). The diffusive resistance is inversely proportional to ðlog nÞ. With each increment in (n), the value of (Rdf) decreases by approximately 37%. The trends for these two parameters are shown in Fig. 5. Both (Rdf(n)) and (Rcn(n)) are non-negative convex functions over [0, 1); therefore, their sum is also a nonnegative convex function with a unique minimum over that range. This minimum can be found via optimization. We now calculate this minimum as well as the value of (n) at which this minimum is achieved.

(15)

We make one more manipulation to this equation, motivated by the fact that the measurement of (dsac) is challenging because it requires histoanatomical analysis of human lactating breast tissue, and it may have a functional relationship to other variables in the equation, in particular, the bifurcation number (n). We make use of a modeling assumption first introduced in Ref. [51] for the bifurcations in the lung, to eliminate (dsac) from this equation. Specifically, if we estimate the maximum space that can be allocated in a lobule with the total length corresponding to infinite number of bifurcations, and furthermore consider the actual space to be represented by the length of the actual bifurcations (n) plus (dsac) that terminates the bifurcations, then one may write dsac þ

n X

Ln ¼

i¼0

1 X

Ln

(16)

i¼0

Substitution for (Ln) in terms of (L0) from Eq. (2) and simple calculations yield

dsac ¼

1 1  2ðnþ1Þ=3 4:85 L0  L ¼ ðnþ1Þ=3 L0 1=3 12 2 1  21=3

071009-4 / Vol. 137, JULY 2015

(17)

Fig. 5 Variation of diffusive and convective resistance versus bifurcation level (subject to average parameter values)

Transactions of the ASME

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

Table 1 The value of various parameters based on direct measurements Parameter Number of lobules per lobe Number of active lobes in each breast Initial duct length Initial duct diameter Blood density Milk density Milk viscosity Blood diffusion coefficient Module of elasticity of alveolar tissue

Symbol

Value

Reference

a b L0 D0 qb qm  Df E

20–40 5–9 8 6 5 (mm) 1.9 6 0.6 (mm) 1060 (kg/m3) 1030 (kg/m3) (1.66–51.90)  106 m2/s (1–12)  1010 (m2/s) 1.95  106 (pa)

[31] [30] [28] [28] [52] [53] [36] [50] [54]

4.1 Optimal Bifurcation Level. For an optimum milk flow in lactating breast, the total resistance which is the summation of the alveolar diffusive resistance in Eq. (18) and the convective resistance in Eq. (11) is minimized. Using the above calculations, the total resistance for the milk flow in a single lobe is RLobe ðcnþdfÞ ¼

128lL0 0:0412E ðn þ 1Þ þ 2n1 pqm D40 2 3  apL0 qb Df

(19)

As mentioned earlier, 90% of breasts have 15–20 lobes; however, only 5–9 lobes in each breast have active mammary ducts that open to the nipple [28]. Since these active ducts are arranged in a parallel configuration, the total resistance in Eq. (19) should be divided by a factor of (b), with 5  b  9. RBreast totðcnþdfÞ ¼

  1 128lL0 2ð2nþ1Þ=3 E ðn þ 1Þ þ b pqm D40 24:25apL0 qb Df

(20)

By taking the derivative with respect to (n) and setting it to zero, it can be verified that the minimum value of the total resistance (Rtot(cn þ ndf)) is achieved as   qm D40 E n ¼ 4:9834 log  18:5721 (21) alL20 qb Df The values of (qm), (qb), (D0), (L0), (), (E), (a), (b), and (Df) are listed in Table 1. These parameters have been determined experimentally, with appropriate references noted in the table. These values depend on the details of the individual anatomy; therefore, they do not have a fixed value across the population and vary between women [28]. Using the average value of each parameter, Eq. (21) leads to the optimum bifurcation level of (n)  25 (Fig. 6).

Recall that the optimization of the bifurcation levels was based on analytical minimization of milk flow resistance. The natural system that is represented by this model is also presumed to be optimized (e.g., by natural selection through generations) but since such a natural optimization is random and may not be mathematically exact, one would expect that our model would yield a value for milk flow resistance that is lower bound on the naturally measured values. Indeed, our model, based on the average measured quantities mentioned above, yields an optimal (minimal) milk flow resistance of logðRtotðcþdfÞ Þ  6:6. To validate our finding, we compared this minimal resistance with the experimental data in the literature. There are several experimental studies measuring the sucking pressure (by infant or breast pump) simultaneously with the milk flow rate during breastfeeding in a healthy breast [33,44,55–61]. The experimentally measured milk flow resistance has been between 7 and 8.1 in log scale (Table 2). The discrepancy of the model-induced value of 6.6 and the experimental range 7–8.1 should be judged in light of the variability of experimental data and the existing assumptions in our mathematical model. For example, assuming the symmetric bifurcation all the way to the alveoli region, expressing the relation between diameter and length of successive generations from Eq. (2), assuming a constant suckling pressure while in the experimental data the sucking pressure has a natural periodic rhythm, neglecting the peripheral pressure exposed to the areola by the infant’s mouth, and possible periductal pressure during suckling. Our model-derived minimum is approximately a factor of 3 below the experimental range; however, the experimental range itself varies by a factor of 12.5 (Fig. 6). Furthermore, the range of experimental milk flow resistance, in the manner calculated in the table, only covers the variation among experiments but not the variation within each experiment. In other words, the resistance for each experiment is reported via an averaging of the pressure and mass production in that experiment. Without such intraexperiment averaging, a slightly higher variation in experimental data will be obtained, as well as a slightly lower gap between the model-based minimum and the experimental data. The empirical reason for this averaging is the periodic nature of infant suckling pressure which produces a broad range in the measurement of

Table 2 Milk flow resistance based on experimental data from the literature

No.

Fig. 6 Optimum bifurcation level (subject to average parameter values)

Journal of Biomechanical Engineering

1 2 3 4 5 6 7 8

Sucking pressure

Mass production

9.5 6 5.3 (cmH2O) 60–160 (cmH2O) 100 (mmHg) 143 6 42 (mmHg) 50 6 5.7 (mmHg) 191.3 6 35.1 (mmHg) 114 6 50 (mmHg) 75 (mmHg)

V_ ¼ 1:28:4 ðml=minÞ V_ ¼ 2836 ðml=minÞ V_ ¼ 5:15  104 ðml=minÞ m_ ¼ 4:6 ðg=5 sÞ m_ ¼ 6:665:9 ðgr=minÞ V_ ¼ 60:663:9 ðml=5 minÞ m_ ¼ 6363:1 ðg=5 minÞ V_ ¼ 43:163:7 ðml=10 minÞ

Log Reference (Rtot) [56] [55] [58] [44] [57] [61] [33] [59]

7.0 7.2 7.3 7.3 7.7 8.0 8.1 8.1

JULY 2015, Vol. 137 / 071009-5

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

Fig. 7 A representation of the Alveolar sacs in the human lactating breast. The immunofluorescence micrograph on the left has been published in part in Hassiotou and Geddes [63]. In the immunohistochemistry micrograph on the right, the two black lines indicate how the two alveolar dimensions were measured.

pressure and flow. One may alternatively report the milk flow resistance measured at the peak as well as the trough of the infant suckling cycle. However, that may produce other methodological issues; therefore, this table resorts to averaging in each experiment. Applying the maximum/minimum values of milk properties and ductal physical parameters (Table 1), results in wide range for 9.2 < dsac < 270.2 lm. To examine how this model-derived range compares with biological data, we measured alveolar diameters in six breast tissues of lactating women obtained from the Archive of the School of Anatomy, Physiology and Human Biology, The University of Western Australia. Paraffinized tissues were sectioned and imaged as described in Ref. [62]. In each of the six tissues, 30 representative alveoli from fuller and emptier adjacent lobules were analyzed. In consistence with the simplification of the presented model of spherical alveoli, complex-shaped alveoli such as the alveolus on the bottom right corner of Fig. 7 (left), were not analyzed. Therefore, the shape of the analyzed alveoli can be described as spherical or oval. For each alveolus, we measured the large and small diameter in a crosslike format (Fig. 7 (right)). The results of the two dimensions for each of the six breast tissues are shown in Fig. 8, which are well within the range of (dsac) determined by our computerized model. By considering the average values of milk properties and ductal physical properties, the total theoretical length of the breast duct can be calculated as

Fig. 8 Variation between and within tissues in alveolar diameter (T1–T6: lactating tissue 1–6). The number of alveoli analyzed per tissue is 30.

071009-6 / Vol. 137, JULY 2015

Ltot ¼

n X

Ln ¼ L0 þ L1 þ …… þ Ln

0

¼ L0 þ ð21=3 ÞL0 þ ð21=3 Þ2 L0 þ …: þ ð21=3 Þn L0 ¼ L0 f1 þ 21=3 þ … þ ð21=3 Þn g ¼ L0

1  ð21=3 Þn  38 mm 1  ð21=3 Þ

(22)

This conclusion is in agreement with and is validated by recent findings based on ultrasound imaging [28] revealing that a large proportion of the glandular tissue appears at 30 mm radius from the base of the nipple. Any separation between the theoretical model and the biological measurements may be explained by the fact that the theoretical model assumes all branches of the bifurcation tree are equally deep, i.e., to (n) levels, while in reality the bifurcation levels in different branches of the duct tree are not always exactly equivalent. Some paths in the actual biological trees are terminated earlier than some others, due to the limitations in the natural amount of space available in the breast for the growth of these trees.

5

Conclusion

Biological organs adapt to their individual tasks in a near-optimal fashion within the constraints of their structure. Our mathematical model uses this principle in the branching structure of human breast to arrive at predictions that are verified by their close match with several clinically measured quantities. The development of this mathematical model has two main aims: (a) understanding of milk flow in the complex geometry of the breast and (b) predictions of the results of experiments which are usually too difficult or costly to carry out directly. Via this model, it is possible to calculate the flow resistance in the breast from the diffusion zone in alveoli (starting point) to removal zone at the nipple openings (end point). This resistance illustrates the relationship between the sucking pressure and milk flow. We also arrive at preliminary formulas that predict parameters such as alveoli dimension and total length of the lobe (from alveolus to the nipple) that vary from subject to subject. Our formulas calculate estimates of these important internal parameters using easily measurable parameters such as milk flow properties and (external) dimensions, e.g., diameter of nipple duct. The prediction accuracy of these formulas has been cross-validated with biopsy samples of alveoli via microscopic imaging. The bifurcation levels in the model are optimized by assuming all the lobules are equally bifurcated. It is found that an optimum structure (i.e., one yielding minimal milk flow resistance) has (n)  25 bifurcations, where the exact value within this range Transactions of the ASME

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

depends on the initial duct size and diameter (in the nipple), as well as other physical properties of the breast. In the present study, minimal needed ingredients are applied to make this model in a way that match with experimental results. Several parameters are chosen as dominant effects (e.g., suckling pressure and milk properties), some as secondary and several items as negligible. Mathematical model is simplified for a greater understanding of the phenomena of milk flow from production source through the nipple in the complex geometry of the breast. Also the resulting model in this paper has the benefit of having fewer unknown parameters, and consequently easier to fit data (fewer experiments are needed). The simplicity of the model is important for two reasons: simpler models are more versatile and also simplicity prevents over-match to the particular data set that is used to derive the parameters of the model. At the same time, the model should admittedly be sophisticated enough so that reasonable fidelity can be attained with respect to the situation at hand. It is reasonable that the initial approach to a model (like the present paper) would produce the simplest form of an eventual family of models. Then, subsequent works can add elaborations to the model and include secondary parameters, and each work can choose the type of secondary parameter that is more suitable for the particular application or a clinical condition that is brought under study. Last but not least, we have focused on a general model for healthy breasts. However, since the formulation is based on the geometry (physical) of the ductal system and milk flow properties (which are measurable), therefore developing a patient specific approach is feasible and beneficial to diagnosis of diseases which are closely related to the ductal system transporting fluids. Future works can build on this work and continue to explore the application of this model in the analysis of various clinical conditions.

Acknowledgment The authors would like to acknowledge Ms. Mary Lee and Ms. Leonie Khoo (School of Anatomy, Physiology and Human Biology, The University of Western Australia) for tissue sectioning; the Translational Pathology Laboratory (TPL) of the University of North Carolina (UNC) at Chapel Hill for expert technical assistance in tissue imaging; and Dr. Ching Tat Lai for assistance in data plotting.

Nomenclature a¼ A¼ b¼ dsac ¼ D¼ Df ¼ E¼ L¼ m_ ¼ DM ¼ n¼ P¼ Dpsac ¼ R¼ u¼ V¼

number of lobules in each lobe cross-sectional area (m) number of active mammary ducts diameter of sac (lm) diameter of duct (mm) diffusion coefficient (m2/s) bulk module of elasticity (Pa) length (mm) mass flow rate (kg/s) deriving potential number of bifurcations pressure (Pa) alveoli muscle contraction pressure (Pa) resistance velocity (m/s) volume (m3)

Greek Symbols d¼ e¼ l¼ ¼ q¼

diffusion path (m) kinetic energy per unit of mass (kJ/kg) dynamic viscosity (Pas) kinematic viscosity (m2/s) density (kg/m3)

Journal of Biomechanical Engineering

Subscript and Superscripts b¼ bf ¼ cn ¼ d¼ df ¼ gen ¼ k at ¼ m¼ n¼ tot ¼

blood bifurcation convection ductal diffusion generation kth bifurcation milk at nth bifurcation total

References [1] Lu, P., Sternlicht, M. D., and Werb, Z., 2006, “Comparative Mechanisms of Branching Morphogenesis in Diverse Systems,” J. Mammary Gland Biol. Neoplasia, 11(3–4), pp. 213–228. [2] Riordan, J., and Wambach, K., 2010, Breastfeeding and Human Lactation, Jones & Bartlett Learning. [3] Love, S. M., and Barsky, S. H., 2004, “Anatomy of the Nipple and Breast Ducts Revisited,” Cancer, 101(9), pp. 1947–1957. [4] Liu, Y., So, R., and Zhang, C., 2002, “Modeling the Bifurcating Flow in a Human Lung Airway,” J. Biomech., 35(4), pp. 465–473. [5] Liu, Y., So, R., and Zhang, C., 2003, “Modeling the Bifurcating Flow in an Asymmetric Human Lung Airway,” J. Biomech., 36(7), pp. 951–959. [6] Yang, X., Liu, Y., and Luo, H., 2006, “Respiratory Flow in Obstructed Airways,” J. Biomech., 39(15), pp. 2743–2751. [7] Kulish, V. V., Lage, J. L., Hsia, C., and Johnson, Jr., R., 2002, “ThreeDimensional, Unsteady Simulation of Alveolar Respiration,” ASME J. Biomech. Eng., 124(5), pp. 609–616. [8] Keyhani, K., Scherer, P., and Mozell, M., 1995, “Numerical Simulation of Airflow in the Human Nasal Cavity,” ASME J. Biomech. Eng., 117(4), pp. 429–441. [9] De Backer, J., Vos, W., Devolder, A., Verhulst, S., Germonpre, P., Wuyts, F., Parizel, P. M., and De Backer, W., 2008, “Computational Fluid Dynamics can Detect Changes in Airway Resistance in Asthmatics After Acute Bronchodilation,” J. Biomech., 41(1), pp. 106–113. [10] Mylavarapu, G., Murugappan, S., Mihaescu, M., Kalra, M., Khosla, S., and Gutmark, E., 2009, “Validation of Computational Fluid Dynamics Methodology Used for Human Upper Airway Flow Simulations,” J. Biomech., 42(10), pp. 1553–1559. [11] Shome, B., Wang, L., Santare, M., Prasad, A., Szeri, A., and Roberts, D., 1998, “Modeling of Airflow in the Pharynx With Application to Sleep Apnea,” ASME J. Biomech. Eng., 120(3), pp. 416–422. [12] Zhao, S., Xu, X., Hughes, A., Thom, S., Stanton, A., Ariff, B., and Long, Q., 2000, “Blood Flow and Vessel Mechanics in a Physiologically Realistic Model of a Human Carotid Arterial Bifurcation,” J. Biomech., 33(8), pp. 975–984. [13] Rayz, V. L., Lawton, M. T., Martin, A. J., Young, W. L., and Saloner, D., 2008, “Numerical Simulation of Pre-and Postsurgical Flow in a Giant Basilar Aneurysm,” ASME J. Biomech. Eng., 130(2), p. 021004. [14] Hughes, P., and How, T., 1996, “Effects of Geometry and Flow Division on Flow Structures in Models of the Distal End-to-Side Anastomosis,” J. Biomech., 29(7), pp. 855–872. [15] Hofer, M., Rappitsch, G., Perktold, K., Trubel, W., and Schima, H., 1996, “Numerical Study of Wall Mechanics and Fluid Dynamics in End-to-Side Anastomoses and Correlation to Intimal Hyperplasia,” J. Biomech., 29(10), pp. 1297–1308. [16] Lei, M., Giddens, D., Jones, S., Loth, F., and Bassiouny, H., 2001, “Pulsatile Flow in an End-to-Side Vascular Graft Model: Comparison of Computations With Experimental Data,” ASME J. Biomech. Eng., 123(1), pp. 80–87. [17] Salzar, R. S., Thubrikar, M. J., and Eppink, R. T., 1995, “Pressure-Induced Mechanical Stress in the Carotid Artery Bifurcation: A Possible Correlation to Atherosclerosis,” J. Biomech., 28(11), pp. 1333–1340. [18] Papaharilaou, Y., Ekaterinaris, J. A., Manousaki, E., and Katsamouris, A. N., 2007, “A Decoupled Fluid Structure Approach for Estimating Wall Stress in Abdominal Aortic Aneurysms,” J. Biomech., 40(2), pp. 367–377. [19] Wang, D., Makaroun, M., Webster, M. W., and Vorp, D. A., 2001, “Mechanical Properties and Microstructure of Intraluminal Thrombus From Abdominal Aortic Aneurysm,” ASME J. Biomech. Eng., 123(6), pp. 536–539. [20] Linninger, A. A., Tsakiris, C., Zhu, D. C., Xenos, M., Roycewicz, P., Danziger, Z., and Penn, R., 2005, “Pulsatile Cerebrospinal Fluid Dynamics in the Human Brain,” IEEE Trans. Biomed. Eng., 52(4), pp. 557–565. [21] Moore, S., David, T., Chase, J., Arnold, J., and Fink, J., 2006, “3d Models of Blood Flow in the Cerebral Vasculature,” J. Biomech., 39(8), pp. 1454–1463. [22] Bonfiglio, A., Leungchavaphongse, K., Repetto, R., and Siggers, J. H., 2010, “Mathematical Modeling of the Circulation in the Liver Lobule,” ASME J. Biomech. Eng., 132(11), p. 111011. [23] Rani, H., Sheu, T. W., Chang, T., and Liang, P., 2006, “Numerical Investigation of Non-Newtonian Microcirculatory Blood Flow in Hepatic Lobule,” J. Biomech., 39(3), pp. 551–563. [24] Debbaut, C., Vierendeels, J., Casteleyn, C., Cornillie, P., Van Loo, D., Simoens, P., Van Hoorebeke, L., Monbaliu, D., and Segers, P., 2012, “Perfusion

JULY 2015, Vol. 137 / 071009-7

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

[25] [26]

[27]

[28]

[29] [30] [31] [32]

[33]

[34]

[35]

[36]

[37] [38] [39] [40] [41] [42] [43] [44]

Characteristics of the Human Hepatic Microcirculation Based on ThreeDimensional Reconstructions and Computational Fluid Dynamic Analysis,” ASME J. Biomech. Eng., 134(1), p. 011003. Cooper, A. P., 1840, On the Anatomy of the Breast, Longman. Going, J. J., and Moffat, D. F., 2004, “Escaping From Flatland: Clinical and Biological Aspects of Human Mammary Duct Anatomy in Three Dimensions,” Pathology, 203(1), pp. 538–544. Ohtake, T., Kimijima, I., Fukushima, T., Yasuda, M., Sekikawa, K., Takenoshita, S., and Abe, R., 2001, “Computer-Assisted Complete Three-Dimensional Reconstruction of the Mammary Ductal/Lobular Systems,” Cancer, 91(12), pp. 2263–2272. Ramsay, D., Kent, J., Hartmann, R., and Hartmann, P., 2005, “Anatomy of the Lactating Human Breast Redefined With Ultrasound Imaging,” J. Anat., 206(6), pp. 525–534. Baum, K. G., McNamara, K., and Helguera, M., 2008, “Design of a Multiple Component Geometric Breast Phantom,” Proc. SPIE, 6913, p. 69134H. Geddes, D. T., 2007, “Inside the Lactating Breast: The Latest Anatomy Research,” J. Midwifery Womens Health, 52(6), pp. 556–563. Sohn, C., Blohmer, J.-U., and Hamper, U. M., 1999, Breast Ultrasound: A Systematic Approach to Technique and Image Interpretation, Thieme. Moffat, D., and Going, J., 1996, “Three Dimensional Anatomy of Complete Duct Systems in Human Breast: Pathological and Developmental Implications,” J. Clin. Pathol., 49(1), pp. 48–52. Geddes, D. T., Kent, J. C., Mitoulas, L. R., and Hartmann, P. E., 2008, “Tongue Movement and Intra-Oral Vacuum in Breastfeeding Infants,” Early Hum. Dev., 84(7), pp. 471–477. Sakalidis, V. S., Kent, J. C., Garbin, C. P., Hepworth, A. R., Hartmann, P. E., and Geddes, D. T., 2013, “Longitudinal Changes in Suck-Swallow-Breathe, Oxygen Saturation, and Heart Rate Patterns in Term Breastfeeding Infants,” J. Hum. Lactation, 29(2), pp. 236–245. Mortazavi, S. N., Geddes, D., and Hassanipour, F., 2014, “Modeling of Milk Flow in Mammary Ducts in Lactating Human Female Breast,” Engineering in Medicine and Biology Society (EMBC), 2014 36th Annual International Conference of the IEEE, pp. 5675–5678. Waller, H., Aschaffenburg, R., and Grant, M. W., 1941, “The Viscosity, Protein Distribution, Andgold Number of the Antenatal and Postnatal Secretions of the Human Mammary Gland,” Biochem. J., 35(3), pp. 272–282. White, F. M., and Corfield, I., 1991, Viscous Fluid Flow, Vol. 2. McGraw-Hill, New York. Bates, J. H., 2009, Lung Mechanics: An Inverse Modeling Approach, Vol. 1. Cambridge University, Cambridge, Wiley, London. Bejan, A., and Lorente, S., 2008, Design With Constructal Theory, Int. J. Engng Ed., 22(1), pp. 140–147, 2006. Bejan, A., 2005, “The Constructal Law of Organization in Nature: Tree-Shaped Flows and Body Size,” J. Exp. Biol., 208(9), pp. 1677–1686. Bejan, A., and Lorente, S., 2006, “Constructal Theory of Generation of Configuration in Nature and Engineering,” J. Appl. Phys., 100(4), p. 041301. Reis, A., Miguel, A., and Aydin, M., 2004, “Constructal Theory of Flow Architecture of the Lungs,” Med. Phys., 31(5), pp. 1135–1140. Bejan, A., and Lorent, S., 2000, Shape and Structure, From Engineering to Nature, Cambridge University, Cambridge, UK. Ramsay, D. T., Mitoulas, L. R., Kent, J. C., Cregan, M. D., Doherty, D. A., Larsson, M., and Hartmann, P. E., 2006, “Milk Flow Rates can be Used to

071009-8 / Vol. 137, JULY 2015

[45]

[46]

[47] [48] [49] [50]

[51]

[52] [53]

[54] [55] [56]

[57]

[58]

[59] [60]

[61]

[62]

[63]

Identify and Investigate Milk Ejection in Women Expressing Breast Milk Using an Electric Breast Pump,” Breastfeed. Med., 1(1), pp. 14–23. Cregan, M. D., Mitoulas, L. R., and Hartmann, P. E., 2002, “Milk Prolactin, Feed Volume and Duration Between Feeds in Women Breastfeeding Their Full-Term Infants Over a 24 h Period,” Exp. Physiol., 87(2), pp. 207–214. Ramsay, D. T., Kent, J. C., Owens, R. A., and Hartmann, P. E., 2004, “Ultrasound Imaging of Milk Ejection in the Breast of Lactating Women,” Pediatrics, 113(2), pp. 361–367. Edgar, A., and Sebring, F., 2005, “Anatomy of a Working Breast,” New Begin., 22(2), pp. 44–50. Neville, M. C., 1990, “The Physiological Basis of Milk Secretiona,” Ann. New York Acad. Sci., 586(1), pp. 1–11. Mills, A. F., 2001, Mass Transfer, Vol. 2. Prentice Hall, Upper Saddle River. Liu, M., Nicholson, J. K., Parkinson, J. A., and Lindon, J. C., 1997, “Measurement of Biomolecular Diffusion Coefficients in Blood Plasma Using Two-Dimensional 1h-1h Diffusion-Edited Total-Correlation NMR Spectroscopy,” Anal. Chem., 69(8), pp. 1504–1509. Nakayama, A., Kuwahara, F., and Sano, Y., 2009, “Why do We Have a Bronchial Tree With 23 Levels of Bifurcation,” J. Heat Mass Transfer, 45(3), pp. 351–354. Comella, M. J., Cutnell, J. D., and Johnson, K. W., 1997, Physics, Wiley. Neville, M. C., Keller, R., Seacat, J., Lutes, V., Neifert, M., Casey, C., Allen, J., and Archer, P., 1988, “Studies in Human Lactation: Milk Volumes in Lactating Women During the Onset of Lactation and Full Lactation,” Am. J. Clin. Nutr., 48(6), pp. 1375–1386. Stacy, R. W., and Giles, F. M., 1959, “Computer Analysis of Arterial Properties,” Circul. Res., 7(6), pp. 1031–1038. Colley, J., and Creamer, B., 1958, “Sucking and Swallowing in Infants,” Br. Med. J., 2(5093), pp. 422–423. Prieto, C., Cardenas, H., Salvatierra, A., Boza, C., Montes, C., and Croxatto, H., 1996, “Sucking Pressure and Its Relationship to Milk Transfer During Breastfeeding in Humans,” J. Reprod. Fertil., 108(1), pp. 69–74. Schrank, W., Al-Sayed, L. E., Beahm, P. H., and Thach, B. T., 1998, “Feeding Responses to Free-Flow Formula in Term and Preterm Infants,” J. Pediatr., 132(3), pp. 426–430. Medoff-Cooper, B., Bilker, W., and Kaplan, J. M., 2010, “Sucking Patterns and Behavioral State in 1-and 2-Day-Old Full-Term Infants,” J. Obstet. Gynecol. Neonat. Nurs., 39(5), pp. 519–524. Drewett, R., and Woolridge, M., 1979, “Sucking Patterns of Human Babies on the Breast,” Early Hum. Dev., 3(4), pp. 315–320. Mitoulas, L. R., Lai, C. T., Gurrin, L. C., Larsson, M., and Hartmann, P. E., 2002, “Effect of Vacuum Profile on Breast Milk Expression Using an Electric Breast Pump,” J. Hum. Lactation, 18(4), pp. 353–360. Mitoulas, L. R., Lai, C. T., Gurrin, L. C., Larsson, M., and Hartmann, P. E., 2002, “Efficacy of Breast Milk Expression Using an Electric Breast Pump,” J. Hum. Lactation, 18(4), pp. 344–352. Hassiotou, F., Hepworth, A. R., Beltran, A. S., Mathews, M. M., Stuebe, A. M., Hartmann, P. E., Filgueira, L., and Blancafort, P., 2013, “Expression of the Pluripotency Transcription Factor oct4 in the Normal and Aberrant Mammary Gland,” Front. Oncol., Vol. 3, Art. 79, pp. 1–15. Hassiotou, F., and Geddes, D., 2013, “Anatomy of the Human Mammary Gland: Current Status of Knowledge,” Clin. Anatomy, 26(1), pp. 29–48.

Transactions of the ASME

Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use