ArticleVolumetric analysis of the terminal ductal lobular unit architecture and cell phenotypes in the human breastGraphical abstractBranch count Proliferation Stromal cells KRT14hi basal epithelium Human breast tissue Volumetric imaging Morphometric analysis Serial sectioning & immunoprofiling Mathematical modelling of branching Highly branched TDLUs are rare in the human breast Consistent shape and general branching parameters Terminal ductal lobular unit (TDLU) Mammary duct TDLU Branching morphogenesis follows evolutionarily conserved principlesHighlightsd Breast terminal ductal lobular unit (TDLU) branching varies regardless of age or parity d TDLUs have amain branch that dominates in bifurcations and has duct-like features d Highly branched TDLUs proliferate more and may contain more S100A4+ stromal cells d Human andmouse mammary glands appear to branch by the same principlesPaavolainen et al., 2024, Cell Reports 43, 114837 October 22, 2024 ª 2024 The Author(s). Published by Elsevier Inc https://doi.org/10.1016/j.celrep.2024.114837Authors Oona Paavolainen, Markus Peurla, Leena M. Koskinen, ..., Colinda L.G.J. Scheele, Pauliina Hartiala, Emilia Peuhu Correspondence emilia.peuhu@utu.fi In brief Paavolainen, Peurla et al. studied intact terminal ductal lobular unit (TDLU) structures to quantitatively assess the branching morphology of the human breast in three dimensions. The study provides a framework for investigating the role of TDLU architecture and cell- type-specific anatomical niches in breast homeostasis and malignant progression.. ll OPEN ACCESS llArticle Volumetric analysis of the terminal ductal lobular unit architecture and cell phenotypes in the human breast Oona Paavolainen,1,2,8 Markus Peurla,1,2,8 Leena M. Koskinen,1,2 Jonna Pohjankukka,1,2 Kamyab Saberi,3 Ella Tammelin,1,2 Suvi-Riitta Sulander,1,2 Masi Valkonen,1 Larissa Mourao,3 Pia Bostro¨m,4,5 Nina Br€uck,4,5 Pekka Ruusuvuori,1 Colinda L.G.J. Scheele,3 Pauliina Hartiala,5,6,7 and Emilia Peuhu1,2,9,* 1Institute of Biomedicine, Cancer Laboratory FICAN West, University of Turku, 20520 Turku, Finland 2Turku Bioscience, University of Turku and A˚bo Akademi University, 20520 Turku, Finland 3VIB Center for Cancer Biology, Department of Oncology, KU Leuven, 3000 Leuven, Belgium 4Department of Pathology, Turku University Hospital, 20520 Turku, Finland 5University of Turku, 20520 Turku, Finland 6Department of Plastic and General Surgery, Turku University Hospital, 20520 Turku, Finland 7Medicity Research Laboratories and InFLAMES Research Flagship Center, University of Turku, 20520 Turku, Finland 8These authors contributed equally 9Lead contact *Correspondence: emilia.peuhu@utu.fi https://doi.org/10.1016/j.celrep.2024.114837SUMMARYThe major lactiferous ducts of the human breast branch out and end at terminal ductal lobular units (TDLUs). Despite their functional and clinical importance, the three-dimensional (3D) architecture of TDLUs has re- mained undetermined. Our quantitative and volumetric imaging of healthy human breast tissue demonstrates that highly branched TDLUs, which exhibit increased proliferation, are uncommon in the resting tissue regardless of donor age, parity, or hormonal contraception. Overall, TDLUs have a consistent shape and branch parameters, and they contain a main subtree that dominates in bifurcation events and exhibits a more duct-like keratin expression pattern. Simulation of TDLU branchingmorphogenesis in three dimensions suggests that evolutionarily conservedmechanisms regulatemammary gland branching in humans andmice despite their anatomical differences. In all, our data provide structural insight into 3D anatomy and branching of the human breast and exemplify the power of volumetric imaging in gaining a deeper understanding of breast biology.INTRODUCTION The human breast epithelium consists of branched ductal net- works formed by 15–20 lobes in which collecting ducts, orig- inating from the nipple, ultimately branch out to form terminal ductal lobular units (TDLUs)1 (Figure 1A). In puberty, systemic hormones, such as growth hormone and estrogen, stimulate the gland to undergo extensive branching morphogenesis, laying down the main anatomical structures of the breast.2 The morphology of the human breast also dramatically changes during the reproductive cycles of pregnancy and lactation in order to adapt to the functional requirements of the tissue. TDLUs are the functional units of the breast, playing a central role during milk production and lactation. Resting TDLUs are typically composed of an extralobular terminal duct (ETD) that branches into lobules with acini at the tips of the ducts (Fig- ure 1A). Interestingly, TDLUs are also considered the predomi- nant anatomical location where breast cancer initiates.4–6 There-Cell Reports 43, 114837, Octo This is an open access article under thefore, a detailed understanding of the structure and composition of the TDLUs is of high clinical importance. While the use of rodent animal models has been pivotal for our basic understanding of healthy mammary gland development and structure, their usefulness for the study of TDLU formation is limited. Rodent mammary glands are formed by a single network of branched ducts within adipose stroma, thus lacking the anatomically distinct lobes and TDLUs within a collagen- rich stroma.7 Whether the branching mechanisms that produce the differential mammary gland outcomes in humans and mice are governed by the same principles has not been previously investigated. Until now, human breast development, anatomy, and three- dimensional (3D) structures have been described mostly based on low-resolution imaging of histological dyes.8,9 While the his- tology and overall anatomy of the human breast epithelium, including the description of different TDLU categories based on donor age and parity, were defined,3 the 3D architecture and exact composition of TDLUs have remained obscure.ber 22, 2024 ª 2024 The Author(s). Published by Elsevier Inc. 1 CC BY license (http://creativecommons.org/licenses/by/4.0/). A B C D E F G H I J (legend on next page) 2 Cell Reports 43, 114837, October 22, 2024 Article ll OPEN ACCESS Article ll OPEN ACCESSAlthough much effort has been put into describing the epithelial and stromal cell populations in the human breast by single-cell analysis,10–13 it is not known how these cell types are anatomi- cally distributed within the breast epithelium or whether their abundance is linked to the extent of TDLU branching. Previous studies assessing cell localization in breast tissue sections7,14,15 or micro-collected TDLUs/ducts16,17 have also provided limited information on the spatial distribution of cell types within TDLUs. In addition, the extent to which TDLU structures vary be- tween individuals and within a breast is unclear. As a result, TDLU composition and structure are still very ill-defined. As recent advances in technology now enable the imaging of larger intact tissue volumes labeled with multiple markers, we employed tissue clearing methods and volumetric light-sheet microscopy to assess the 3D morphology of TDLUs derived from healthy breast tissue donors, both parous and nulliparous, in different age categories. We quantified TDLU structures and related their size and complexity to age, parity, and hormonal contraception status. Volumetric imaging was followed by serial sectioning and multiplex immunofluorescence analysis of the same TDLUs to relate cell type distribution and abundance to their 3D structure. Finally, we combined our quantitative analysis with mathematical modeling to investigate whether general branching patterns or principles exist for the human TDLU. RESULTS Volumetric imaging of optically cleared human breast tissue reveals that highly branched TDLUs are rare regardless of donor age or parity The morphology of entire human TDLUs has previously been documented only with bright-field or confocal microscopy imag- ing, thus lacking detailed volumetric information of the struc- ture.3,18 In the early studies conducted on large numbers of breast tissue wholemounts using low-resolution imaging, TDLUs were divided into four different branching types, where lobular type 1 represents young nulliparous tissue (6–17 acini/ TDLU), type 2 mature nulliparous tissue (17–60 acini/TDLU), type 3 parous tissue (>65 acini/TDLU), and type 4 branching occurring only during pregnancy and lactation.3 More recently, however, the number of acini per TDLU was reported to be com- parable between nulliparous cases and samples collected 18 months19 or 5 years20 postpartum.Figure 1. Volumetric imaging of optically cleared human breast tissue or parity (A) A schematic representation of the human breast terminal ductal lobular unit (T human breast TDLU (donor 5). Scale bar: 100 mm. (B) A schematic of the sample preparation and mounting protocol. (C) The age, parity, and use of hormonal contraceptives by the tissue donors. Th (D) Data analysis pipeline. A light-sheet microscopy image of a TDLU visualized w color coding (middle), and an image of the branching structure model generated structure starting point. Scale bars: 100 mm. Determination of branch count per (E) Representative images of breast tissue after iFLASH clearing and light-sheet (F) The number of branches per TDLU with respect to previous lobular typing.3 L (G) The number of branches per TDLU graphed in relation to donor age. Parity of t 15 donors). (H–J) The number of branches per TDLU graphed in relation to donor age and pa Whitney U test), or days since last menstruation started (J) (n = 26 TDLUs from 6To assess the influence of age, parity, and hormonal contra- ception on the TDLU branch count using volumetric imaging, we obtained breast tissue from healthy premenopausal women of different age (18–48) and parity (0–4 children) undergoing breast reduction surgery (Figures 1B and 1C; Table S1). Tissue pieces of approximately 5–20 mm in diameter were pre- screened for the presence of mammary epithelium by nuclear la- beling, cleared with the FUnGI (fructose, urea, and glycerol for imaging) protocol18 or our improved version of the FLASH (fast light-microscopic analysis of antibody-stained whole organs) optical tissue clearing protocol21 (termed hereon iFLASH), and labeled by immunofluorescence for cytokeratin 8 (KRT8) to visu- alize luminal mammary epithelial cells (Figure 1B). We obtained good sample transparency and were able to visualize entire TDLUs (Figure S1; Video S1). Altogether, 67 TDLUs from 15 do- nors were imaged by volumetric light-sheet imaging (Figure 1C) (deconvoluted raw data are available at the BioImage Archive: https://doi.org/10.6019/S-BIAD1136). From the 3D image stacks, TDLU structures were analyzed to visualize and quantify the branched networks (Figure 1D). The KRT8+ luminal epithe- lium formed a continuous tubular structure that bifurcated in a seemingly random manner and ended as enlarged acini (Figures 1D and 1E), demonstrating that the imaging protocol could reveal intact TDLU structures. Quantification of the total 3D branch count for all analyzed TDLUs indicated that most TDLUs were medium sized, corre- sponding roughly to type 2 lobules of an earlier categorization3 (Figure 1F). The sample cohort was strongly divided into young (18–28) nulliparous and mature (35–48) parous (>6 years post- partum or unknown) donors, likely due to the relatively high average age of first-time mothers in Finland (30.1 years)22 (Figures 1C and 1G; Table S1). Importantly, the abundance of highly branched TDLUs was not higher in the samples from mature parous donors (Figures 1G and 1H). In turn, the heteroge- neity in TDLU branch count was the largest in samples from young nulliparous donors (Figure 1G). More than half of the do- nors were using hormonal contraceptives, but no emergent ef- fect on TDLU branch count could be attributed to hormonal contraception (Figure 1I; Table S1). Although based on a limited number of samples from women not using hormonal contracep- tion, a higher TDLU branch count was detected in the breasts of women in themid-phase of their menstrual cycle (Figure 1J). This is in agreement with previous studies demonstrating higherreveals that highly branched TDLUs are rare regardless of donor age DLU) based on current knowledge and an H&E-stained histological section of e number of TDLUs analyzed per donor is indicated. ith depth color coding (KRT8, left), a segmented image visualized with the same with Imaris FilamentTracer (right). Cyan dot in the image indicates the branch TDLU demonstrated with illustration. microscopy. Scale bars: 100 mm. ob1: 0–17 branches, Lob2: 35–60 branches, and Lob3: >65 branches. he donor or the use of hormonal contraceptives is indicated (n = 67 TDLUs from rity (H), hormonal contraceptive use (I) (n = 67 TDLUs from 15 donors; Mann- donors; mean ± SEM with non-linear curve fit). Cell Reports 43, 114837, October 22, 2024 3 A C D E F G H B Figure 2. TDLUs contain a main subtree that dominates in bifurcation events and exhibits a more duct-like KRT14 expression pattern (A and B) Confocal microscopy imaging of KRT14 (basal epithelial, cyan), KRT8 (luminal epithelial, red), and nuclei (blue) in immunolabeled frozen human breast tissue sections (A) or in entire TDLUs of cleared breast tissue pieces (B). Higher levels of KRT14 are detected in the extralobular terminal duct (ETD) compared to the lobules of TDLUs. Regions of interest (ROIs) depict the ETDs (ROI 1) and lobuli (ROI 2). Scale bars: (A) 100 mm (main) and 25 mm (ROIs); (B) 250 mm (main) and 50 mm (ROIs). n = 6–8 from 3 donors. (legend continued on next page) 4 Cell Reports 43, 114837, October 22, 2024 Article ll OPEN ACCESS Article ll OPEN ACCESSestrogen levels, and thus elevated mammary epithelial prolifera- tion, approximately at this cycle phase.23,24 In all, we generated a workflow for clearing and volumetric imaging of human breast tissue. Our data suggest that parity, in association with higher age, is not linked to a higher number of TDLU branches (or acini) in the human breast several years postpartum, thus confirming previous reports demonstrating similar findings based on tissue sections.19,20 TDLUs contain a main subtree that dominates in bifurcation events and exhibits a more duct-like KRT14 expression pattern The location of different cell types in human breast tissue has been investigated using tissue sections7,14,15 or in micro- collected TDLUs/ducts.16,17 However, the cellular composition or tissue-level organization in different parts of the TDLU remain poorly understood. KRT8 and cytokeratin 14 (KRT14) intermedi- ate filament proteins have been established as luminal and pre- dominantly basal epithelial markers, respectively, in the human mammary gland.13,25 Immunolabeling of KRT8 and KRT14 in healthy breast tissue sections revealed higher KRT14 expression in the basal cell layer of the ETDs compared to the lobules (Fig- ure 2A), as previously observed.15,26 In order to build a more comprehensive view of the cellular or- ganization of TDLUs, we cleared breast tissue pieces with iFLASH protocol and immunolabeled for KRT8 and KRT14. The samples were imaged by confocal (Figure 2B) or light-sheet mi- croscopy imaging (Figures 2C and 2D; Video S2). Although some KRT14 expression was detected in the entire basal mammary epithelial cell layer (Figures 2C and 2D), a subpopulation of basal cells expressing high levels of KRT14 (KRT14hi) was detected predominantly in the ETDs of whole TDLUs (Figures 2C and 2D), as expected based on the tissue section analysis (Figure 2A). Interestingly, volumetric imaging and quantitative analysis of the KRT14hi basal cells present in the individual subtrees (Figures 2E–2H) revealed that the amount of KRT14hi cells re- duces gradually from the ETD toward the acini. Furthermore, the main subtree of a TDLU, defined as the subtree with the most bifurcation events (levels), presents with a higher number of KRT14hi basal cells per branch level (Figures 2G and S2). In other words, the main subtree of each TDLU was significantly more likely to harbor the highest number of KRT14hi cells per branch level (Figure 2H). Conversely, the expression of the luminal marker KRT8 seemingly increased in intensity toward the acini (Figures 2B–2D). Immunolabeling of another basal marker protein, transgelin, produced a clearly different expres- sion pattern than KRT14, where expression levels were highest at the lobular level around the acini and lower close to the ETD (Figures S3A and S3B). In all, these data reveal an additional level(C and D) Light-sheet microscopy imaging of KRT14 (basal epithelial, cyan) and K (C). KRT14hi cells are visible in the duct (arrow), ETDs (ROI 1), and lobuli (ROI 2). Lo Scale bars: 100 mm (main) and 50 mm (ROIs). Images are representative of 28 TD (E and F) Determination of subtrees in TDLU branching structures (E). Represent branching hierarchy (TDLU 5, cyan; line width is representative of the number of (G) Quantification of the branch levels and KRT14hi cells in individual subtrees of (n = 28, from 9 donors). Separate subtrees defined to begin at branch level = 3 a (H) Quantification of the presence of KRT14hi cells in main vs. not-main subtreesof organization within the TDLU structure, identifying amain sub- tree that dominates in bifurcation events and exhibits a more duct-like KRT14 expression pattern when compared to the smaller subtrees of the same TDLU. Multiplex IHC of serially sectioned TDLUs demonstrates higher cell proliferation in highly branched TDLUs In addition to the heterogeneity in KRT patterns, we observed variation in TDLU size and branching complexity, with branch counts ranging from 8 to 208 per TDLU (Figures 1F and S2). To provide insight into the currently unknown tissue-level fea- tures associated with different degrees of branching in human TDLUs (Figure 1), we conducted serial sectioning of 110 TDLUs from samples previously imaged by light-sheet micro- scopy (Figures 1 and 2) and performed multiplex immunohisto- chemistry (IHC) using antibodies for known mammary gland cell markers (Figure 3A). H&E staining was performed at 40 mm intervals (every 10th tissue section) for the entire tissue sample volume (Figure 3B), which could also be processed into a 3D reconstruction histology image (Figure 3C). Calculation of the cumulative branch count from all the H&E sections of each TDLU (Figures 3D and S4) resulted in a comparable distribution of the branch counts per TDLU to what was detected by light- sheet imaging (Figure 1F), even though the number of TDLUs that could be visualized from the serial sections was higher. Cat- egories from low to high TDLU branching were assigned based on the quartiles of branch count per TDLU in the whole dataset (1st–4th quartiles are equivalent to ranges of 0%–25%, 25%– 50% 50%–75%, and 75%–100% in branch count, respectively) (Figure 3D). The OPAL 6-plex immunolabeling panels containing estab- lished markers for the mammary epithelium (E-cadherin [CDH1], pan-keratin), luminal epithelial cells (KRT8, cytokeratin 19 [KRT19], estrogen receptor [ER]), basal myoepithelial cells (p63, a-smooth muscle actin [aSMA, data not shown]),15,27 T cells (CD3), stromal cells (S100A4, also known as FSP-1),28 and cell proliferation (Ki67)29 were labeled at 2–3 z-positions per sample array (Figure 3E). The samples were imaged and in- dividual TDLUs analyzed quantitatively for the relative amount of selected markers in the epithelial or stromal areas. The level of each marker was correlated with the parity and age of the donor (Figure 3F) and the quartiles of TDLU branch count (Figure 3G). ER was expressed in a subset of KRT19+ luminal cells, as ex- pected,30 and exhibited no clear correlation with the TDLU branch count, although a trend toward lower numbers of ER+ luminal cells was observed in highly branched TDLUs (Figures 3E and 3G). In turn, the amount of ER+ cells was higher in the mammary epithelium of mature parous donors compared to young nulliparous (Figure 3F), as was previously shown toRT8 (luminal epithelial, red) in entire TDLUs within cleared breast tissue pieces wer level of KRT14 is detected in the rest of the basal mammary epithelium (D). LUs from 9 donors. ative tree structures of TDLUs with the distribution of KRT14hi cells within the KRT14hi cells in each branch; branch length not to scale) (F). the TDLUs demonstrating that the dominant subtrees retain high KRT14 levels nd persist for a minimum of 3 levels after cutoff. (n = 27–31 subtrees; Fisher’s exact test). Cell Reports 43, 114837, October 22, 2024 5 A C B E D F G H J K 1st quartile 4th quartile KRT19 ER KRT19 ER PanCK Ki67 PanCK Ki67 KRT8 CD3 KRT8 CD3 CDH1 S100A4 CDH1 S100A4 KRT19 p63 KRT19 p63 Labelling Serial sectioning Serial H&E staining (step size = 40μm) Tissue microarrays Light sheet samples Labelling Labelling Histology Histology Labelling ER Ki67 p63 CD3 (intra-TDLU stroma) CD3 (epithelium) S100A4 (intra-TDLU stroma) S100A4 (epithelium) 0.0000 0.0005 0.0010 0.0015 1st 2nd 3rd 4th 0.0000 0.0005 0.0010 0.0015 0.0020 1st 2nd 3rd 4th 0.00 0.02 0.04 0.06 0.08 1st 2nd 3rd 4th 1st 2nd 3rd 4th 1st 2nd 3rd 4th 0.00 0.05 0.10 0.15 0.20 1st 2nd 3rd 4th 0.000 0.005 0.010 0.015 0.020 0.025 1st 2nd 3rd 4th R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l R el at iv e le ve l 0.000 0.002 0.004 0.006 0.0008 0.1482 0.0251 NP <30 P >35 0.000 0.001 0.002 0.003 0.004 ns NP <30 P >35 NP <30 P >35 NP <30 P >35 NP <30 NP P >35 P NP P 0.0 0.1 0.2 0.3 0.4 ns 0.00 0.05 0.10 0.15 ns 0.0 0.1 0.2 0.3 0.4 0.5 ns ns 0 100 200 300 B ra nc h co un t 4th 3rd 2nd 1st quartile 0.00 0.02 0.04 0.06 0.08 0.10 0.0315 0.00 0.02 0.04 0.06 0.08 0.00 0.02 0.04 0.06 0.08 0.000 0.005 0.010 0.015 0.020 0.025 μm μm 320 160 240 250 250 500 500 80 0 0 0 μm I (legend on next page) 6 Cell Reports 43, 114837, October 22, 2024 Article ll OPEN ACCESS Article ll OPEN ACCESSoccur during aging.31 On the contrary, p63+ basal myoepithelial cells were detected in equal amounts in the human TDLUs irre- spective of age, parity, or TDLU branch count (Figures 3E–3G). Interestingly, proliferation (Ki67+ cells) was higher in TDLUs with the highest branch count (4th quartile) compared to TDLUs with a lower branch count (2nd quartile) (Figures 3E–3G), sug- gesting that highly branched TDLUs could be in a more prolifer- ative state compared to an average-sized TDLU. In the stromal compartment of the TDLU, we observed CD3+ T cells in equal amounts irrespective of age, parity, or TDLU branch count (Figures 3E–3G). Interestingly, mature and parous women had more intraepithelial CD3+ T cells than young nullip- arous donors (Figure 3H), but the amount of intraepithelial CD3+ T cells was not correlated to TDLUbranch count (Figure 3I). An interesting trend, although not statistically significant, was observed in stromal S100A4+ cells that include fibroblasts, T cells, monocytes, and macrophages.28 S100A4 labeling was somewhat higher in the stroma of TDLUs with the highest branch count (4th quartile) compared to TDLUs with a lower branch count (Figures 3E–3G). Few S100A4+ cells were detected within the epithelium regardless of TDLU branch count (Figures 3J and 3K). The finding that highly branched TDLUs appear more popu- lated by S100A4+ cells is well in line with previous studies demonstrating that the presence of fibroblasts and macro- phages promotes the branching morphogenesis of mammary epithelial cells in vitro32,33 and in vivo in mice.34 The exact nature of the S100A4+ cells remains to be investigated further. Together, our data provide evidence of the clinical relevance of stromal cells’ regulatory role in the branching morphogenesis of the human TDLU. Volumetric imaging of optically cleared human breast tissue reveals consistent branching parameters between different TDLUs and individual donors The relatively high variation in TDLU branch counts prompted us to investigate the morphological similarity of individual branches in TDLUs. In order to quantitatively establish the 3D morphology and dimensions of the human TDLU,weperformed tissue clearing with the FUnGI protocol that, unlike several other protocols, has been reported to preserve the original size of the specimen.18,35 Fourteen KRT8-labeled TDLU structures from four donors were successfully imaged with volumetric light-sheet microscopy (Fig- ure 4A). The donors represented different ages (range: 24–50 years) and parity (range: 0–2 children) (Table S1; Figure S5A). For quantitative branching analysis, the 3D image stacks were segmented for KRT8 signals, followed by quantitative analysis of the TDLU structures (Figures 4B and 4C; Video S1).Figure 3. Multiplex immunohistochemistry of serially sectioned TDLUs (A) Serial section sample preparation protocol. Tissue samples utilized for volume serially sectioned (4-mm sections). Every 10th section was H&E stained, allowing f used for multiplex labeling. (B–D) An example of H&E series of a TDLU (B). 3D reconstruction histology of the s on branch count calculated from the H&E-stained serial sections (n = 110 from 1 (E) OPAL 6-plex immunolabeling with antibodies against ER (cyan), p63 (red), Ki6 25 mm (ROIs). (F–K) Quantification of specific IHC signals/cell types per TDLU in function of dono 10 donors; mean ± SEM; Mann-Whitney U test).The morphometric analysis of human breast TDLUs demon- strates regularity in branch length, diameter, and angle (Figures 4D, 4E, and S5C–S5E). The average branch length was 71 mm (SD: 37.8 mm), the average diameter was 35 mm (SD: 8.9 mm), and the average branch angle with respect to the parent duct was 49 (SD: 29.9) (Figure 4E). The TDLU tip branches (acini) were, on average, 1.2 times as wide as the ducts in the TDLU (Figure 4F). A TDLU formed most commonly from 4–6 consecutive bifurcation events, and branches at levels higher than 10 consecutive bifurcations were few (Figure 4G). Average TDLU length (distance from ETD to the farthest acinus; Figure 4H) and width (largest width in the direction perpendic- ular to the ETD direction; Figure 4H) were 460 and 360 mm, respectively, and TDLUs typically had a regular shape with an aspect ratio of 1.5 (Figure 4H). These measurements are similar to the previously reported TDLU sizes in H&E tissue sections.36 Interestingly, the branching of benign mammary epithelium adjacent to a preinvasive lesion (ductal carcinoma in situ) was comparable to healthy donor tissues (Figure 4D, donor 16). In all, our data demonstrate that TDLUs have a consistent shape and that their branch parameters are largely comparable be- tween different TDLUs and individual donors of variable age and parity (Figure 4D). In contrast to previous observations us- ing tissue sections,37 the total number of branches or acini (tips) per TDLUwas not proportional to the TDLU diameter (Figure 4I), suggesting that branching density may vary. Human mammary gland branching morphogenesis appears to be governed by the same principles as in mouse Our volumetric analysis highlighted that despite several com- mon features, such as the existence of a dominating TDLU sub- tree (Figure 2), and the consistent TDLU shape and branch parameters (Figure 4), breast TDLU branching is irregular, asymmetric, and variable in branch count (Figures S5F, S5G, and S6). To investigate the potential mechanisms of TDLU for- mation, we considered existing models for branching morpho- genesis. A fractal branching model with symmetrically occur- ring bifurcations and gradually decreasing branch length and diameter was previously proposed for human TDLUs,38 but we found this model to be incompatible with the 3D branching structure of the TDLU (Figures S5C–S5E). The fact that human TDLU branching networks are unique and vary greatly in size (Figures S4, S5F, and S5G) also speaks against a stereotyped branching morphogenesis as observed in the lung39 or kid- ney.40 Additionally, the absence of loop structures does not support models where the final tree-like form is reacheddemonstrates higher cell proliferation in highly branched TDLUs tric light-sheet imaging were re-embedded into tissue microarrays (TMAs) and or 3D reconstruction of tissue anatomy. Remaining, consecutive sections were tructure in (C). ETD (orange arrow). (D) Distribution of all analyzed TDLUs based 0 donors). Scale bars: (B) 200 mm. 7 (magenta), CD3 (orange), and S100A4 (green). Scale bars: 50 mm (main) and r age/parity (F, H, and J) or percentiles of branch count (G, I, and K) (n = 110 from Cell Reports 43, 114837, October 22, 2024 7 A B D C 200 400 600 100 200 300 400 500 TDLU length (μm) TD LU w id th (μ m ) E F G di am et er 1 2 3 16 1 2 3 16 1 2 3 16 H I Average TDLU diameter (μm) Donor Donor Donor ROI Figure 4. Volumetric imaging of optically cleared human breast tissue reveals consis- tent branching parameters between different TDLUs and individual donors (A and B) Representative image of breast tissue cleared with FUnGI, imaged (KRT8) with light- sheet microscopy (A) and the branch architecture segmented for quantitative branch analysis (B). Scale bar: 100 mm. (C) Parameter definitions. Top image: TDLU length (red line) and TDLU width (green line), cyan dot in- dicates the network starting point (ETD). Bottom image: branch length (red line), branch diameter (green line), and branching angle (blue). (D) Distributions of the parameters branch length, branch diameter, and branching angle by donor. Mean values of individual TDLUs are shown by the white dots, and the gray area shows value distri- bution. (E) Pooled distributions of branch length, diameter, and angle (nlength = 689, ndiameter = 689, and nangle = 689 branches from 14 TDLUs of 4 donors). (F) Distributions of acinar and ductal diameters per TDLU (n = 14 from 4 donors). (G) Distribution of branching levels. Data are pooled from all complete TDLU structures for which the whole network could be reliably traced (n = 8 TDLUs from 4 donors). (H) Measured length and diameter of the whole TDLU. Line shows linear fit to data, excluding the outlier point. TDLU aspect ratio from the fit is 1.47 (n = 10 TDLUs from 4 donors). (I) Total number of branches and tip branches (acini) in TDLU as a function of average TDLU diameter. Average TDLU diameter is the average of the measured TDLU length and width (n = 10 TDLUs from 4 donors). In (E) and (F), the boxes show median and interquartile range, and the whiskers show the full range of data points. Article ll OPEN ACCESSthrough remodeling of a mesh of interconnected ducts, like in the pancreas.41 Branching morphogenesis by branching and annihilating random walk (BARW) has previously been used to simulate the generation of mouse mammary gland,42,43 kidney,42 and pancreas,44 zebrafish nervous system,45 and cholangiocyte orga- noids.46 Whether human TDLUs branch and develop with the same principles as the mouse mammary gland has not been ad- dressed. In the BARW model, branches grow forward, bifurcate randomly, and terminate if they meet an existing branch.42 Growing branches try to avoid collisions with other branches.45 To study the human TDLU, we adapted the original 2D BARW model42 to three dimensions and utilized our experimentally ob- tained parameters for branch diameter, bifurcation probability, termination distance, repulsion distance, and branching and rota- tion angle distributions (Figures 4E, 5A, S5H, and S5I; Table S2). In the simulations, growth takes placemostly at the perimeter of the network by branches invading empty space, and the networks will grow indefinitely unless growth is somehow limited.We tested two mechanisms for limiting network growth: (1) growth was limited to a predefined volume42 or (2) with a model in which branches have gradually waning growth potential that eventually stops the branching. The latter model was inspired by the accu- mulation of KRT14hi cells into one TDLU subtree that reached higher branch levels (Figure 2), which could mirror the existence of unequal and potentially reducing growth capacity.8 Cell Reports 43, 114837, October 22, 2024In the volume-limited model, the network was allowed to grow within a volume equal to the experimentally observed average TDLU volume. Any branch reaching the edges of the simulation volume was terminated. We used an adaptive limiting volume (Figure 5B) that follows the network growth direction tomaximize the use of available space by preventing occasional premature termination by all branches colliding directly with network vol- ume edges. The volume-limited model produced ductal struc- tures very similar to the TDLUs in human tissue samples (Figures 5E–5J and S6). Simulated networks had similar branch counts to those observed experimentally (Figure 5F), and the networks were equally effectively filling the space (Figure 5D). Model prediction of branch-level distribution followed the gen- eral trend of that observed experimentally (R2 = 0.99, Figure 5G). The model also predicts well the fraction of terminal branches (acini) of the total number of branches (Figure S5B). The topology, ancestry of the subtrees, and connectedness of the networks were compared as before.42,46 The probability of a branch terminating at a given level (Figure 5H) follows a similar trend observed before,42,46 leveling to approximately 0.6 at levels above 5. Leveling happens at earlier levels compared to mousemammary glands due to the lower branch count in human TDLUs. Simulated networks follow the trend well (R2 = 0.90, Fig- ure 5H). Tree plots of the branch ancestry (Figure S6) show that human TDLU subtrees are not self-similar or symmetric, but instead their size varies substantially within a TDLU. Many BBranching free zone radius Termination distance Repulsion distance Branch diameter Volume limited model Gradual waning model SimulationExperimentExperiment Simulation A B C D E F G H I J K L M N O P Figure 5. Simulation of human TDLU branching morphogenesis with volume-limited BARW model resembles more accurately the experi- mental data (A) Schematic image drawn to scale visualizing simulation parameters. (B) Simulation volume control in volume-limited model. Red and green dots mark the network origin and center of mass, respectively. Dashed line indicates the general direction of the network, defined by a vector from network origin to center of mass. (C) Simulation volume control in gradual waning model. Main image shows the fraction of KRT14hi cells on a level of all the KRT14hi cells in the TDLU. Dots mark the average value and error bars the variation of all TDLUs. Inset shows the fit of the exponential equation to the experimental KRT14hi cell distribution used in the model as the waning boundary between branching (orange area) and self-termination (blue area). n = 28 TDLUs of 9 donors. (legend continued on next page) Cell Reports 43, 114837, October 22, 2024 9 Article ll OPEN ACCESS Article ll OPEN ACCESSshallow subtrees are typically observed, and the highest branch- ing occurs in 1–2 ancestral subtree lines that continue growing up to levels 10 and above (Figure S6). Comparable behavior can also be generated in the volume-limited model (Figures S6, 5I, and 5J). The fit of the model results to experimental observa- tions is good for subtree persistence (R2 = 0.99, Figure 5I) and subtree size in branches (R2 = 0.99, Figure 5J). Overall the 3D simulation based on the BARW model with limiting volume pro- duces a close fit to the experimental data. In the gradually waning branching potential model, branch behavior gradually shifts from branching to self-termination as level increases. Self-termination represents a situation where the branch ceases to grow, no more branching is instigated, and the branch becomes inactive. Hypothesizing that KRT14hi cells could represent a feature of more duct-like cells with higher branching potential, we estimated the probability boundary be- tween branching and self-termination from the experimentally observed decline of the KRT14hi cell distribution in the TDLUs (Figure 5C). Themodelofgraduallywaningbranchingpotentialwas lesssuc- cessful in reproducing the features of the experimental networks. In this model, the shape of the networks was more spherical or even flatter, with an aspect ratio <1 (Figure 5K). We introduced and optimized another model parameter, the global guidance force,45 to reproduce the realistic network shapewithout affecting the topology results. In general, this model also produced a satis- factory fit to the experimental data (Figures 5L–5P). However, these networks grew much larger in volume, as indicated by the much less densely packed structure (Figure 5D). The model was not able to produce networks in which subtrees would persist to levels above 15 (Figure 5O) and subtrees consisted of fewer branches (Figure 5P) even if the overall branch count distribution was similar to that observed experimentally (Figure 5L). Growth in these networks takes place in several shallow subtrees growing in parallel rather than in 1–2 main subtrees that extend to high branching levels, as was observed experimentally, and is better modeled by the volume limiting model (Figures 5 and S6). In summary, the simulations demonstrate that branching morphogenesis of human TDLUs appears to be governed by the same principles as in mouse, with the volume-limited 3D BARW model capturing well the morphology of the experimen- tally analyzed TDLU structures. DISCUSSION The 3D branching structure of the human breast has remained largely unexplored, particularly in high resolution and in a quan-(D) Branching density for experimental TDLUs and both simulationmodels. Branch and it is given in units 1/(106 mm3). nexperiment = 8 TDLUs of 4 donors and nsimulati (E and K) Images of experimental and simulated networks for both models. Blue (F and L) Simulation predictions for the number of branches in a TDLU compared (G and M) Distribution of branching levels in experimental data and simulation re (H and N) Probability of branch termination at different branching levels. (I and O) Probability of a subtree beginning at level 3 to span a number of levels. (J and P) Probability of a subtree beginning at level 3 to contain a number of bra In (G)–(J) and (M)–(P), dots and error bars indicate mean and standard deviation standard deviation of results of 1,000 simulations. For (F)–(J) and (L)–(P), nexperim In (D), (F), and (L), the boxes show median and interquartile range, and the whisk 10 Cell Reports 43, 114837, October 22, 2024titative manner. Here, we set out to provide structural insight into the 3D anatomy, cell type location, and branching of resting hu- man TDLUs. TDLUs of different donors were heterogeneous in branch count but had a consistent shape and branch length, width, and angle. Our data provide evidence that resting TDLUs of young nulliparous women are indistinguishable from TDLUs of mature parous donors in terms of branch count but have lower amounts of ER+ luminal epithelial cells and intraepi- thelial T cells. While hormonal contraception did not seem to in- fluence the overall branch count of TDLUs in our sample cohort, branching varied during the menstrual cycle, as expected. Highly branched TDLUs had more proliferating epithelial cells compared to less branched TDLUs and, interestingly, exhibited a trend toward a higher amount of stromal cells surrounding the epithelium. We simulated branchingmorphogenesis of the human TDLUby adapting the BARW model42 to three dimensions and using experimentally defined parameters for the modeling. Two poten- tial mechanisms were considered for limiting the network growth: (1) limitation of growth to a predefined volume, as previously used in the modeling of mouse mammary gland branching,42 or (2) lim- itation of growth by gradually waning growth potential. The latter mechanism was based on the gradual decline observed in KRT14hi cells toward the acini of the TDLU. KRT14 has previously been linked to increased adhesive and migratory properties,47,48 as well as the mechanical integrity of cells.49 Thus, we hypothe- sized that KRT14 expression level could potentially reflect local differences in growth or migration potential and presumed that a decline in branching potential could provide a sufficient mecha- nism to regulate TDLU growth. While the simulated networks re- sulting from mathematical modeling with this mechanism for limiting the network growth had similarities to the experimental data, they were sparse and lacked a dominating subtree. Instead, the networks simulated with the volume-limited BARW model highly resembled the TDLUs in the clinical spec- imen despite the fairly wide variation in the data, suggesting that human mammary gland branching morphogenesis is gov- erned by the same principles as the mouse mammary gland. The most significant difference between the species in terms of the modeling parameters is the branching area, which in hu- mans was limited to a single TDLU. Our data point to the fact that the TDLU form is likely limited by interaction with the stroma. Why TDLUs evolved and why they adopt this particular size and shape in humans remain unknown. Furthermore, whether the morphogenesis of an entire lobe would follow comparable prin- ciples with TDLU branching remains to be studied. Specifically, it would be interesting to address if the larger branch diameter ining density is calculated as the number of branches per whole network volume, on = 1,000 for each model. dot indicates the network origin. Line and point sizes are not to scale. to experimental results. sults. nches. of the experimental data, and the orange line and blue shaded area mean and ent = 63 TDLUs of 15 donors and nsimulation = 1,000 for each model. ers show the full range of data points. Article ll OPEN ACCESSbigger mammary ducts is sufficient to account for the higher branch length. Terminatedmammary ductal tips typically reside close to an ex- isting duct or the fat pad boundary in mice.50 Additionally, ductal branching was shown to occur in vitro only when the tip was remote from the other ducts.51 Consequently, ductal elongation and branching was suggested to proceed as a ‘‘default state,’’ with tip termination occurring only when the tips come into prox- imity to existing ducts or the edges of the fat pad.42 This scenario, as well as the modeling presented here, suggests that growth limiting mechanisms exist around the TDLUs as well. While there is currently limited information on the negative regulation of branching morphogenesis, transforming growth factor b1 (TGF- b1) secretion51,52 and progressive envelopment of the gland by extracellular matrix53,54 have been suggested to limit mammary gland branching. Therefore, gradual alterations in TDLU-proximal stroma may contribute to TDLU growth inhibition. Mammary gland stromal cells have been shown to promote mammary epithelial branching.32–34We observed a trend toward an elevated amount of S100A4+ stromal cells in the TDLUs of the highest branch count where proliferation was also higher than in less branched TDLUs. S100A4 expression has been detected, for instance, in fibroblasts55–58 and macrophages.59,60 Two distinct fibroblast lineages, lobular and interlobular, have been described in the human breast.61 Intralobular fibroblasts, also quantified here, have been shown to support epithelial growth and morphogenesis.61 Fibroblast growth factor signaling was previously demonstrated to regulate the ability of stromal fibroblasts to induce mouse mammary epithelial branching in vitro.32 In addition, perturbed hedgehog signaling62 or integrin activity regulation63 in the mouse mammary gland stroma was shown to disturb mammary gland branching morphogenesis. In light of these data, it is plausible that elevated amounts of in- tralobular fibroblasts in highly branched, more proliferative TDLUs could be functionally linked to the elevated branch count. However, such causality remains to be confirmed. In earlier studies conducted with lower-resolution imaging, TDLUs were divided into four different types of branching.3 Of the TDLU samples analyzed quantitatively in this study, only a small fraction contained more than 81 branches in total, thus reaching the category Lob3 that was earlier suggested to predominate in parous women.3 Importantly, our data suggest that the TDLUs of young nulliparous and mature parous do- nors are, in fact, not significantly different in branch counts. Similar conclusions have also been made by others who have investigated human TDLUs with tissue sections and concluded that the length of the postpartum period is an important factor affecting branching.19,20 In our sample cohort, the shortest known postpartum period to sample collection was 6 years. Although the dataset analyzed here is still relatively small and lacks the ability to distinguish the ef- fect of age from parity, our data suggest that the effect of par- ity on TDLU branching may not be as straightforward as pre- viously described. In addition to parity and age, menstrual cycling64 and hormon- al contraception3 have been indicated to affect TDLU branching. While no differences in the TDLU branch count could be clearly linked to the hormonal contraception status in our study, anelevated branch count was detected in the samples from donors in mid-cycle, corroborating earlier findings demonstrating increased proliferation at a similar phase.64–66 This increase is followed by a decrease toward the end of the cycle, which matches the reported increase in apoptosis.67,68 Further, the po- sition of the TDLU within the whole gland could affect the branching type, as fibroglandular density, along with tumor inci- dence, has been indicated to be highest at the upper outer quad- rant of the breast in Western women.69 The samples in our branch count analysis were collected from reduction mammo- plasty operations, where the tissue is resected from the lower quadrants of the breast. Higher ER expression was observed in the TDLU epithelium of samples from mature parous women compared to young nullip- arous women. Gradually increasing ER status has been previ- ously reported in response to aging,30,31 and ER+ breast cancer ismore commonwith increasing age.70 In turn, parity and the use of hormonal contraception have been shown to decrease ER levels.71–73 In our sample cohort, where the effects of age and parity cannot be separated, the effect of aging appears to domi- nate over parity. Also, T cell infiltration into the mammary epithe- lium was more prominent in TDLUs of mature parous donors. Similar effects of age or parity have not been previously reported, and a more detailed characterization of this T cell population is warranted. Overall, this study presents a systematic analysis and mathe- matical modeling of the glandular 3D structure in healthy human breast epithelium. Our findings add an additional layer of detail to breast anatomy and demonstrate that 3D histology or volumetric imaging can be used to study the relationship between TDLU ar- chitecture and the cell types that form it. This study also provides a reference dataset for future development of light-sheet image analysis and a protocol for future investigations on the role of ductal architecture in the regulation of growth in healthy andma- lignant breast epithelium, including carcinogenesis and invasive progression. Importantly, the identities and niches of distinct cell types within the TDLU structure are likely to become clearer as research proceeds, revealing structure-function relationships in the breast tissue.Limitations of the study As the volumetric imaging and immunoprofiling approaches used in this study are currently not compatible with the analysis of very large sample sets, the data are based on a limited cohort of 15 healthy tissue donors. The donors were either young nulliparous or mature parous, thus limiting any conclu- sions related to the individual contributions of these factors. As such, our findings are an important addition to the current understanding of human breast anatomy since thorough 3D morphometric analysis of human breast anatomy has not been conducted before.RESOURCE AVAILABILITY Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Emilia Peuhu (emilia. peuhu@utu.fi).Cell Reports 43, 114837, October 22, 2024 11 Article ll OPEN ACCESSMaterials availability This study did not generate new unique reagents. Data and code availability Raw light-sheet imaging data were deposited on the BioImage Archive: https://doi.org/10.6019/S-BIAD1136. The mathematical model used for the 3D simulation of TDLU branching is available on Mendeley Data: https://doi. org/10.17632/y8d52ygdxv.1.74 Any additional information required to reana- lyze the data reported in this paper is available from the lead contact upon request. ACKNOWLEDGMENTS We thank the patients that voluntarily donated their tissue to this study; the clinical staff of Turku University Hospital and Mediclinic, Oud-Heverlee (BE), for assisting in sample collection; Minna Santanen, Ciaran Butler-Hallissey, Barbara Ramos Artigot, and Jukka Karhu for technical assistance; and all the lab members for constructive feedback over the course of this project. The Cell Imaging and Cytometry Core facility (Turku Bioscience, University of Turku, A˚bo Akademi University, and Biocenter Finland), Biomedicum Imag- ing Unit (University of Helsinki), Histocore (Institute of Biomedicine, University of Turku), and Turku Bioimaging are acknowledged for services, instrumenta- tion, and expert advice. This work was supported by the Sigrid Juse´lius Foun- dation, an Academy of Finland research fellowship (323096), and the Hospital District of Southwest Finland (11083). O.P. has been supported by the Turku Doctoral Program ofMolecular Medicine. L.M. would like to acknowledge sup- port from Stichting Tegen Kanker (fellowship number ZKE4757). AUTHOR CONTRIBUTIONS Conceptualization, O.P., M.P., and E.P.; methodology, O.P., M.P., M.V., L.M., P.R., C.L.G.J.S., and E.P.; formal analysis, O.P., M.P., L.M.K., J.P., and M.V.; investigation, O.P., M.P., L.M.K., J.P., K.S., E.T., S.-R.S., M.V., L.M., N.B., P.B., and P.H.; writing – original draft, O.P., M.P., and E.P.; writing – review & editing, all authors; visualization, O.P., M.P., and M.V.; supervision, P.R., C.L.G.J.S., P.H., and E.P.; funding acquisition, C.L.G.J.S. and E.P. DECLARATION OF INTERESTS The authors declare no competing interests. STAR+METHODS Detailed methods are provided in the online version of this paper and include the following: d KEY RESOURCES TABLE d EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILSB Human tissue samples d METHOD DETAILS B Tissue clearing B Histology and immunohistochemistry B Light sheet microscopy and deconvolution B Confocal imaging B Image segmentation and branching analysis B Paraffin embedding of cleared tissue samples B 3D reconstruction B Multiplex staining B Multiplex image analysis B Mathematical modeling d QUANTIFICATION AND STATISTICAL ANALYSIS SUPPLEMENTAL INFORMATION Supplemental information can be found online at https://doi.org/10.1016/j. celrep.2024.114837.12 Cell Reports 43, 114837, October 22, 2024Received: March 23, 2024 Revised: August 20, 2024 Accepted: September 20, 2024 REFERENCES 1. Bannister, L., Berry, M., Collins, P., Dyson, M., Dussek, J., and Ferguson, M. (1995). Gray’s Anatomy: The anatomical basis of medicine and surgery. In Gray’s Anatomy (Churchill), pp. 417–424. 2. Macias, H., and Hinck, L. (2012). Mammary gland development. WIREs Developmental Biology 1, 533–557. https://doi.org/10.1002/wdev.35. 3. Russo, J., Rivera, R., and Russo, I.H. (1992). Influence of age and parity on the development of the human breast. Breast Cancer Res. Treat. 23, 211–218. https://doi.org/10.1007/BF01833517. 4. Wellings, S.R., Jensen, H.M., and Marcum, R.G. (1975). An Atlas of Sub- gross Pathology of the Human Breast With Special Reference to Possible Precancerous Lesions2. J. Natl. Cancer Inst. 55, 231–273. https://doi.org/ 10.1093/jnci/55.2.231. 5. Russo, J., and Russo, I.H. (1997). Differentiation and breast cancer. Me- dicina 57, 81–91. 6. Taba´r, L., Dean, P.B., Tucker, F.L., Yen, A.M.-F., Fann, J.C.-Y., Lin, A.T.-Y., Smith, R.A., Duffy, S.W., and Chen, T.H.-H. (2022). Breast cancers originating from the terminal ductal lobular units: In situ and invasive acinar adenocarcinoma of the breast. AAB. European Journal of Radiology 152, 110323. https://doi.org/10.1016/j.ejrad.2022.110323. 7. Dontu, G., and Ince, T.A. (2015). Of Mice and Women: A Comparative Tis- sue Biology Perspective of Breast Stem Cells and Differentiation. J. Mammary Gland Biol. Neoplasia 20, 51–62. https://doi.org/10.1007/ s10911-015-9341-4. 8. Russo, J., Romero, A.L., and Russo, I.H. (1994). Architectural pattern of the normal and cancerous breast under the influence of parity. Cancer Ep- idemiol. Biomarkers Prev. 3, 219–224. 9. Howard, B.A., and Gusterson, B.A. (2000). Human Breast Development. J. Mammary Gland Biol. Neoplasia 5, 119–137. https://doi.org/10.1023/ A:1026487120779. 10. Nguyen, Q.H., Pervolarakis, N., Blake, K., Ma, D., Davis, R.T., James, N., Phung, A.T., Willey, E., Kumar, R., Jabart, E., et al. (2018). Profiling human breast epithelial cells using single cell RNA sequencing identifies cell diversity. Nat. Commun. 9, 2028. https://doi.org/10.1038/s41467-018- 04334-1. 11. Pal, B., Chen, Y., Vaillant, F., Capaldo, B.D., Joyce, R., Song, X., Bryant, V.L., Penington, J.S., Di Stefano, L., Tubau Ribera, N., et al. (2021). A sin- gle-cell RNA expression atlas of normal, preneoplastic and tumorigenic states in the human breast. EMBO J. 40, e107333. https://doi.org/10. 15252/embj.2020107333. 12. Chen, Y., Pal, B., Lindeman, G.J., Visvader, J.E., and Smyth, G.K. (2022). R code and downstream analysis objects for the scRNA-seq atlas of normal and tumorigenic human breast tissue. Sci. Data 9, 96. https://doi.org/10. 1038/s41597-022-01236-2. 13. Gray, G.K., Li, C.M.-C., Rosenbluth, J.M., Selfors, L.M., Girnius, N., Lin, J.-R., Schackmann, R.C.J., Goh, W.L., Moore, K., Shapiro, H.K., et al. (2022). A human breast atlas integrating single-cell proteomics and tran- scriptomics. Dev. Cell 57, 1400–1420.e7. https://doi.org/10.1016/j.dev- cel.2022.05.003. 14. Villadsen, R., Fridriksdottir, A.J., Rønnov-Jessen, L., Gudjonsson, T., Rank, F., LaBarge, M.A., Bissell, M.J., and Petersen, O.W. (2007). Evi- dence for a stem cell hierarchy in the adult human breast. JCB (J. Cell Biol.) 177, 87–101. https://doi.org/10.1083/jcb.200611114. 15. Kumar, T., Nee, K., Wei, R., He, S., Nguyen, Q.H., Bai, S., Blake, K., Pein, M., Gong, Y., Sei, E., et al. (2023). A spatially resolved single-cell genomic atlas of the adult human breast. Nature 620, 181–191. https://doi.org/10. 1038/s41586-023-06252-9. Article ll OPEN ACCESS16. Goldhammer, N., Kim, J., Villadsen, R., Rønnov-Jessen, L., and Petersen, O.W. (2022). Myoepithelial progenitors as founder cells of hyperplastic hu- man breast lesions upon PIK3CA transformation. Commun. Biol. 5, 219. https://doi.org/10.1038/s42003-022-03161-x. 17. Kohler, K.T., Goldhammer, N., Demharter, S., Pfisterer, U., Khodosevich, K., Rønnov-Jessen, L., Petersen, O.W., Villadsen, R., and Kim, J. (2022). Ductal keratin 15+ luminal progenitors in normal breast exhibit a basal- like breast cancer transcriptomic signature. npj Breast Cancer 8, 81. https://doi.org/10.1038/s41523-022-00444-8. 18. Rios, A.C., Capaldo, B.D., Vaillant, F., Pal, B., van Ineveld, R., Dawson, C.A., Chen, Y., Nolan, E., and Fu, N.Y.; 3DTCLSM Group (2019). Intraclo- nal Plasticity in Mammary Tumors Revealed through Large-Scale Single- Cell Resolution 3D Imaging. Cancer Cell 35, 618–632.e6. https://doi.org/ 10.1016/j.ccell.2019.02.010. 19. Jindal, S., Gao, D., Bell, P., Albrektsen, G., Edgerton, S.M., Ambrosone, C.B., Thor, A.D., Borges, V.F., and Schedin, P. (2014). Postpartum breast involution reveals regression of secretory lobules mediated by tissue-re- modeling. Breast Cancer Res. 16, R31. https://doi.org/10.1186/bcr3633. 20. Ogony, J., de Bel, T., Radisky, D.C., Kachergus, J., Thompson, E.A., Deg- nim, A.C., Ruddy, K.J., Hilton, T., Stallings-Mann, M., Vachon, C., et al. (2022). Towards defining morphologic parameters of normal parous and nulliparous breast tissues by artificial intelligence. Breast Cancer Res. 24, 45. https://doi.org/10.1186/s13058-022-01541-z. 21. Messal, H.A., Almagro, J., Zaw Thin, M., Tedeschi, A., Ciccarelli, A., Blackie, L., Anderson, K.I., Miguel-Aliaga, I., van Rheenen, J., and Beh- rens, A. (2021). Antigen retrieval and clearing for whole-organ immunoflu- orescence by FLASH. Nat. Protoc. 16, 239–262. https://doi.org/10.1038/ s41596-020-00414-z. 22. Perinatal statistics - parturients, delivers and newborns - THL (2024). Finnish Institute for Health and Welfare (THL), Finland. https://thl.fi/en/ statistics-and-data/statistics-by-topic/sexual-and-reproductive-health/ parturients-deliveries-and-births/perinatal-statistics-parturients-delivers- and-newborns. 23. Potten, C.S., Watson, R.J., Williams, G.T., Tickle, S., Roberts, S.A., Harris, M., and Howell, A. (1988). The effect of age and menstrual cycle upon pro- liferative activity of the normal human breast. Br. J. Cancer 58, 163–170. https://doi.org/10.1038/bjc.1988.185. 24. Laidlaw, I.J., Clarke, R.B., Howell, A., Owen, A.W., Potten, C.S., and An- derson, E. (1995). The proliferation of normal human breast tissue im- planted into athymic nude mice is stimulated by estrogen but not proges- terone. Endocrinology 136, 164–171. https://doi.org/10.1210/endo.136.1. 7828527. 25. Taylor-Papadimitriou, J., Stampfer, M., Bartek, J., Lewis, A., Boshell, M., Lane, E.B., and Leigh, I.M. (1989). Keratin expression in human mammary epithelial cells cultured from normal andmalignant tissue: relation to in vivo phenotypes and influence of medium. J. Cell Sci. 94, 403–413. https://doi. org/10.1242/jcs.94.3.403. 26. Nguyen, Q.H., Pervolarakis, N., Blake, K., Ma, D., Davis, R.T., James, N., Phung, A.T., Willey, E., Kumar, R., Jabart, E., et al. (2018). Profiling human breast epithelial cells using single cell RNA sequencing identifies cell di- versity. Nat. Commun. 9, 2028. https://doi.org/10.1038/s41467-018- 04334-1. 27. Prabhu, S.D., Rai, H.S., Nayak, R., Naik, R., and Jayasheelan, S. (2023). Study of the Immunohistochemical Expression of p63 in Benign Lesions and Carcinoma of the Breast at a Tertiary Hospital in South India. Cureus 15, e48557. https://doi.org/10.7759/cureus.48557. 28. Fendt, B.M., Hirschmann, A., Bruns, M., Camarillo-Retamosa, E., Ospelt, C., and Vogetseder, A. (2023). Protein atlas of fibroblast specific protein 1 (FSP1)/S100A4. Histol. Histopathol. 38, 1391–1401. https://doi.org/10. 14670/HH-18-621. 29. van Dierendonck, J.H., Keijzer, R., van de Velde, C.J., and Cornelisse, C.J. (1989). Nuclear distribution of the Ki-67 antigen during the cell cycle: com- parison with growth fraction in human breast cancer cells. Cancer Res. 49, 2999–3006.30. Santandrea, G., Bellarosa, C., Gibertoni, D., Cucchi, M.C., Sanchez, A.M., Franceschini, G., Masetti, R., and Foschini, M.P. (2021). Hormone Recep- tor Expression Variations in Normal Breast Tissue: Preliminary Results of a Prospective Observational Study. J. Personalized Med. 11, 387. https:// doi.org/10.3390/jpm11050387. 31. Shoker, B.S., Jarvis, C., Sibson, D.R., Walker, C., and Sloane, J.P. (1999). Oestrogen receptor expression in the normal and pre-cancerous breast. J. Pathol. 188, 237–244. https://doi.org/10.1002/(SICI)1096-9896(199907) 188:3<237::AID-PATH343>3.0.CO;2-8. 32. Sumbal, J., and Koledova, Z. (2019). FGF signaling in mammary gland fi- broblasts regulates multiple fibroblast functions and mammary epithelial morphogenesis. Development 146, dev185306. https://doi.org/10.1242/ dev.185306. 33. Sumbal, J., Fre, S., and Sumbalova Koledova, Z. (2024). Fibroblast- induced mammary epithelial branching depends on fibroblast contrac- tility. PLoS Biol. 22, e3002093. https://doi.org/10.1371/journal.pbio. 3002093. 34. Gouon-Evans, V., Rothenberg, M.E., and Pollard, J.W. (2000). Postnatal mammary gland development requires macrophages and eosinophils. Development 127, 2269–2282. 35. Almagro, J., Messal, H.A., Zaw Thin, M., van Rheenen, J., and Behrens, A. (2021). Tissue clearing to examine tumour complexity in three dimen- sions. Nat. Rev. Cancer 21, 718–730. https://doi.org/10.1038/s41568- 021-00382-w. 36. Figueroa, J.D., Pfeiffer, R.M., Patel, D.A., Linville, L., Brinton, L.A., Gierach, G.L., Yang, X.R., Papathomas, D., Visscher, D., Mies, C., et al. (2014). Ter- minal Duct Lobular Unit Involution of the Normal Breast: Implications for Breast Cancer Etiology. J. Natl. Cancer Inst. 106, dju286. https://doi. org/10.1093/jnci/dju286. 37. Rosebrock, A., Caban, J.J., Figueroa, J., Gierach, G., Linville, L., Hewitt, S., and Sherman, M. (2013). Quantitative Analysis of TDLUs using Adaptive Morphological Shape Techniques. Proc. SPIE 8676, 86760N. https://doi.org/10.1117/12.2006619. 38. Honeth, G., Schiavinotto, T., Vaggi, F., Marlow, R., Kanno, T., Shinomiya, I., Lombardi, S., Buchupalli, B., Graham, R., Gazinska, P., et al. (2015). Models of Breast Morphogenesis Based on Localization of Stem Cells in the Developing Mammary Lobule. Stem Cell Rep. 4, 699–711. https:// doi.org/10.1016/j.stemcr.2015.02.013. 39. Metzger, R.J., Klein, O.D., Martin, G.R., and Krasnow, M.A. (2008). The branching programme of mouse lung development. Nature 453, 745–750. https://doi.org/10.1038/nature07005. 40. Lefevre, J.G., Short, K.M., Lamberton, T.O., Michos, O., Graf, D., Smyth, I.M., and Hamilton, N.A. (2017). Branching morphogenesis in the devel- oping kidney is governed by rules that pattern the ureteric tree. Develop- ment 144, 4377–4385. https://doi.org/10.1242/dev.153874. 41. Dahl-Jensen, S.B., Yennek, S., Flasse, L., Larsen, H.L., Sever, D., Karre- more, G., Novak, I., Sneppen, K., and Grapin-Botton, A. (2018). Decon- structing the principles of ductal network formation in the pancreas. PLoS Biol. 16, e2002842. https://doi.org/10.1371/journal.pbio.2002842. 42. Hannezo, E., Scheele, C.L.G.J., Moad, M., Drogo, N., Heer, R., Sampo- gna, R.V., van Rheenen, J., and Simons, B.D. (2017). A Unifying Theory of Branching Morphogenesis. Cell 171, 242–255.e27. https://doi.org/10. 1016/j.cell.2017.08.026. 43. Nerger, B.A., Jaslove, J.M., Elashal, H.E., Mao, S., Kosmrlj, A., Link, A.J., and Nelson, C.M. (2021). Local accumulation of extracellular matrix regu- lates global morphogenetic patterning in the developing mammary gland. Curr. Biol. 31, 1903–1917.e6. https://doi.org/10.1016/j.cub.2021.02.015. 44. Sznurkowska, M.K., Hannezo, E., Azzarelli, R., Rulands, S., Nestorowa, S., Hindley, C.J., Nichols, J., Go¨ttgens, B., Huch,M., Philpott, A., and Simons, B.D. (2018). Defining Lineage Potential and Fate Behavior of Precursors during Pancreas Development. Dev. Cell 46, 360–375.e5. https://doi. org/10.1016/j.devcel.2018.06.028.Cell Reports 43, 114837, October 22, 2024 13 Article ll OPEN ACCESS45. Uc¸ar, M.C., Kamenev, D., Sunadome, K., Fachet, D., Lallemend, F., Ada- meyko, I., Hadjab, S., and Hannezo, E. (2021). Theory of branching morphogenesis by local interactions and global guidance. Nat. Commun. 12, 6830. https://doi.org/10.1038/s41467-021-27135-5. 46. Roos, F.J.M., van Tienderen, G.S., Wu, H., Bordeu, I., Vinke, D., Albarinos, L.M., Monfils, K., Niesten, S., Smits, R., Willemse, J., et al. (2022). Human branching cholangiocyte organoids recapitulate functional bile duct for- mation. Cell Stem Cell 29, 776–794.e13. https://doi.org/10.1016/j.stem. 2022.04.011. 47. Fujiwara, S., Deguchi, S., andMagin, T.M. (2020). Disease-associated ker- atin mutations reduce traction forces and compromise adhesion and col- lective migration. J. Cell Sci. 133, jcs243956. https://doi.org/10.1242/jcs. 243956. 48. Cheung, K.J., Gabrielson, E., Werb, Z., and Ewald, A.J. (2013). Collective Invasion in Breast Cancer Requires a Conserved Basal Epithelial Program. Cell 155, 1639–1651. https://doi.org/10.1016/j.cell.2013.11.029. 49. Ramms, L., Fabris, G., Windoffer, R., Schwarz, N., Springer, R., Zhou, C., Lazar, J., Stiefel, S., Hersch, N., Schnakenberg, U., et al. (2013). Keratins as themain component for themechanical integrity of keratinocytes. Proc. Natl. Acad. Sci. USA 110, 18513–18518. https://doi.org/10.1073/pnas. 1313491110. 50. Silberstein, G.B. (2001). Postnatal mammary gland morphogenesis. Microsc. Res. Tech. 52, 155–162. https://doi.org/10.1002/1097- 0029(20010115)52:2<155::AID-JEMT1001>3.0.CO;2-P. 51. Nelson, C.M., VanDuijn, M.M., Inman, J.L., Fletcher, D.A., and Bissell, M.J. (2006). Tissue Geometry Determines Sites of Mammary Branching Morphogenesis in Organotypic Cultures. Science 314, 298–300. https:// doi.org/10.1126/science.1131000. 52. Silberstein, G.B., and Daniel, C.W. (1987). Reversible Inhibition of Mam- mary Gland Growth by Transforming Growth Factor-b. Science 237, 291–293. https://doi.org/10.1126/science.3474783. 53. Silberstein, G.B., and Daniel, C.W. (1984). Glycosaminoglycans in the basal lamina and extracellular matrix of serially aged mouse mammary ducts. Mech. Ageing Dev. 24, 151–162. https://doi.org/10.1016/0047- 6374(84)90067-8. 54. Osin, P.P., Anbazhagan, R., Bartkova, J., Nathan, B., and Gusterson, B.A. (1998). Breast development gives insights into breast disease. Histopa- thology 33, 275–283. https://doi.org/10.1046/j.1365-2559.1998.00479.x. 55. Cheng, N., Bhowmick, N.A., Chytil, A., Gorksa, A.E., Brown, K.A., Mur- aoka, R., Arteaga, C.L., Neilson, E.G., Hayward, S.W., and Moses, H.L. (2005). Loss of TGF-b type II receptor in fibroblasts promotes mammary carcinoma growth and invasion through upregulation of TGF-a-MSP- and HGF-mediated signaling networks. Oncogene 24, 5053–5068. https://doi.org/10.1038/sj.onc.1208685. 56. Trimboli, A.J., Cantemir-Stone, C.Z., Li, F., Wallace, J.A., Merchant, A., Creasap, N., Thompson, J.C., Caserta, E., Wang, H., Chong, J.-L., et al. (2009). Pten in Stromal Fibroblasts Suppresses Mammary Epithelial Tu- mors. Nature 461, 1084–1091. https://doi.org/10.1038/nature08486. 57. Pickup, M.W., Hover, L.D., Polikowsky, E.R., Chytil, A., Gorska, A.E., Nov- itskiy, S.V., Moses, H.L., and Owens, P. (2015). BMPR2 loss in fibroblasts promotes mammary carcinoma metastasis via increased inflammation. Mol. Oncol. 9, 179–191. https://doi.org/10.1016/j.molonc.2014.08.004. 58. Koledova, Z., Zhang, X., Streuli, C., Clarke, R.B., Klein, O.D., Werb, Z., and Lu, P. (2016). SPRY1 regulates mammary epithelial morphogenesis by modulating EGFR-dependent stromal paracrine signaling and ECM re- modeling. Proc. Natl. Acad. Sci. USA 113, E5731–E5740. https://doi. org/10.1073/pnas.1611532113. 59. Le Hir, M., Hegyi, I., Cueni-Loffing, D., Loffing, J., and Kaissling, B. (2005). Characterization of renal interstitial fibroblast-specific protein 1/S100A4- positive cells in healthy and inflamed rodent kidneys. Histochem. Cell Biol. 123, 335–346. https://doi.org/10.1007/s00418-005-0788-z. 60. O¨sterreicher, C.H., Penz-O¨sterreicher, M., Grivennikov, S.I., Guma, M., Koltsova, E.K., Datz, C., Sasik, R., Hardiman, G., Karin, M., and Brenner,14 Cell Reports 43, 114837, October 22, 2024D.A. (2011). Fibroblast-specific protein 1 identifies an inflammatory sub- population of macrophages in the liver. Proc. Natl. Acad. Sci. USA 108, 308–313. https://doi.org/10.1073/pnas.1017547108. 61. Morsing, M., Klitgaard, M.C., Jafari, A., Villadsen, R., Kassem, M., Pe- tersen, O.W., andRønnov-Jessen, L. (2016). Evidence of two distinct func- tionally specialized fibroblast lineages in breast stroma. Breast Cancer Res. 18, 108. https://doi.org/10.1186/s13058-016-0769-2. 62. Zhao, C., Cai, S., Shin, K., Lim, A., Kalisky, T., Lu, W.-J., Clarke, M.F., and Beachy, P.A. (2017). Stromal Gli2 activity coordinates a niche signaling program for mammary epithelial stem cells. Science 356, eaal3485. https://doi.org/10.1126/science.aal3485. 63. Peuhu, E., Kaukonen, R., Lerche, M., Saari, M., Guzma´n, C., Rantakari, P., De Franceschi, N., Wa¨rri, A., Georgiadou, M., Jacquemet, G., et al. (2017). SHARPIN regulates collagen architecture and ductal outgrowth in the developing mouse mammary gland. EMBO J. 36, 165–182. https://doi. org/10.15252/embj.201694387. 64. Ramakrishnan, R., Khan, S.A., and Badve, S. (2002). Morphological Changes in Breast Tissue with Menstrual Cycle. Mod. Pathol. 15, 1348– 1356. https://doi.org/10.1097/01.MP.0000039566.20817.46. 65. Going, J.J., Anderson, T.J., Battersby, S., andMacIntyre, C.C. (1988). Pro- liferative and secretory activity in human breast during natural and artificial menstrual cycles. Am. J. Pathol. 130, 193–204. 66. So¨derqvist, G., Isaksson, E., Schoultz, B. von, Carlstro¨m, K., Tani, E., and Skoog, L. (1997). Proliferation of breast epithelial cells in healthy women during the menstrual cycle. Am. J. Obstet. Gynecol. 176, 123–128. https://doi.org/10.1016/S0002-9378(97)80024-5. 67. Anderson, T.J., Ferguson, D.J., and Raab, G.M. (1982). Cell turnover in the ‘‘resting’’ human breast: influence of parity, contraceptive pill, age and lat- erality. Br. J. Cancer 46, 376–382. 68. Longacre, T.A., and Bartow, S.A. (1986). A correlative morphologic study of human breast and endometrium in the menstrual cycle. Am. J. Surg. Pathol. 10, 382–393. https://doi.org/10.1097/00000478- 198606000-00003. 69. Chen, J.-H., Liao, F., Zhang, Y., Li, Y., Chang, C.-J., Chou, C.-P., Yang, T.-L., and Su, M.-Y. (2017). 3D MRI for Quantitative Analysis of Quadrant Percent Breast Density: Correlation with Quadrant Location of Breast Cancer. Acad. Radiol. 24, 811–817. https://doi.org/10.1016/j.acra.2016. 12.016. 70. Shah, A., Haider, G., Abro, N., Bhutto, S., Baqai, T.I., Akhtar, S., and Ab- bas, K. (2022). Correlation Between Age and Hormone Receptor Status in Women With Breast Cancer. Cureus 14, e21652. https://doi.org/10. 7759/cureus.21652. 71. Russo, J., Ao, X., Grill, C., and Russo, I.H. (1999). Pattern of distribution of cells positive for estrogen receptor a and progesterone receptor in relation to proliferating cells in the mammary gland. Breast Cancer Res. Treat. 53, 217–227. https://doi.org/10.1023/A:1006186719322. 72. Williams, G., Anderson, E., Howell, A., Watson, R., Coyne, J., Roberts, S.A., and Potten, C.S. (1991). Oral contraceptive (OCP) use increases pro- liferation and decreases oestrogen receptor content of epithelial cells in the normal human breast. Int. J. Cancer 48, 206–210. https://doi.org/10. 1002/ijc.2910480209. 73. Hilton, H.N., Clarke, C.L., and Graham, J.D. (2018). Estrogen and proges- terone signalling in the normal breast and its implications for cancer devel- opment. Mol. Cell. Endocrinol. 466, 2–14. https://doi.org/10.1016/j.mce. 2017.08.011. 74. Peurla, M. (2023). 3D branching simulation v. 3.21. 1. https://doi.org/10. 17632/y8d52ygdxv.1. 75. Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., et al. (2012). Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682. https://doi.org/10.1038/nmeth.2019. 76. Arganda-Carreras, I., Kaynig, V., Rueden, C., Eliceiri, K.W., Schindelin, J., Cardona, A., and Sebastian Seung, H. (2017). Trainable Weka Article ll OPEN ACCESSSegmentation: a machine learning tool for microscopy pixel classification. Bioinformatics 33, 2424–2426. https://doi.org/10.1093/bioinformatics/ btx180. 77. Short, K., Hodson, M., and Smyth, I. (2013). Spatial mapping and quanti- fication of developmental branching morphogenesis. Development 140, 471–478. https://doi.org/10.1242/dev.088500. 78. Bankhead, P., Loughrey, M.B., Ferna´ndez, J.A., Dombrowski, Y., McArt, D.G., Dunne, P.D., McQuaid, S., Gray, R.T., Murray, L.J., Coleman, H.G., et al. (2017). QuPath: Open source software for digital pathology image anal- ysis. Sci. Rep. 7, 16878. https://doi.org/10.1038/s41598-017-17204-5. 79. Woodford, O. (2024). vol3d v2 - File Exchange - MATLAB Central. https:// se.mathworks.com/matlabcentral/fileexchange/22940-vol3d-v2.80. Susaki, E.A., Tainaka, K., Perrin, D., Kishino, F., Tawara, T., Watanabe, T.M., Yokoyama, C., Onoe, H., Eguchi, M., Yamaguchi, S., et al. (2014). Whole-Brain Imaging with Single-Cell Resolution Using Chemical Cock- tails and Computational Analysis. Cell 157, 726–739. https://doi.org/10. 1016/j.cell.2014.03.042. 81. Ruusuvuori, P., Valkonen, M., Kartasalo, K., Valkonen, M., Visakorpi, T., Nykter, M., and Latonen, L. (2022). Spatial analysis of histology in 3D: quantification and visualization of organ and tumor level tissue envi- ronment. Heliyon 8, e08762. https://doi.org/10.1016/j.heliyon.2022. e08762. 82. Vince, J. (2021). Vector Analysis for Computer Graphics (Springer). https:// doi.org/10.1007/978-1-4471-7505-6.Cell Reports 43, 114837, October 22, 2024 15 Article ll OPEN ACCESSSTAR+METHODSKEY RESOURCES TABLEREAGENT or RESOURCE SOURCE IDENTIFIER Antibodies Rat monoclonal anti-keratin 8 (KRT8) (clone TROMA-I) antibody (Tissue pieces; 1:100, IHC-F; 1:1000) DSHB Cat#TROMA-I-c; RRID:AB_531826 Rabbit polyclonal anti-keratin 14 (KRT14) (clone Poly19053) antibody (Tissue pieces; 1:100, IHC-F; 1:1000) Biolegend Cat#905301; RRID:AB_2565048 Rabbit polyclonal anti-transgelin (TAGLN) antibody (Tissue pieces; 1:50, IHC-F; 1:200) Abcam Cat#ab155272; RRID:AB_2637003 Mouse monoclonal anti-p63 (clone D-9) antibody (mIHC; 1:500) Santa Cruz Cat#sc-25268; RRID:AB_628092 Rabbit monoclonal anti-estrogen receptor a (clone D6R2W) antibody (mIHC; 1:50) Cell Signaling Cat#13258S; RRID:AB_2632959 Mouse monoclonal anti-keratin 19 (KRT19) (clone A53-B/A2.26) antibody (mIHC; 1:1000) Thermo Fisher Scientific Cat#MA5-12663; RRID:AB_10984317 Rabbit monoclonal anti-FSP1/S100A4 antibody (mIHC; 1:1200) Sigma-Aldrich Cat#SAB5701199; RRID: AB_3532201 Rabbit monoclonal anti-E Cadherin (phospho S838 + S840) (clone [EP913(2)Y]) antibody (mIHC; 1:3000) Abcam Cat#ab76319; RRID:AB_2076796 Rat monoclonal anti-keratin 8 (KRT8) (clone TROMA-I) antibody (mIHC; 1:1000) Sigma-Aldrich Cat#MABT329; RRID:AB_2891089 Rat monoclonal anti-CD3E (clone CD3-12) antibody (mIHC; 1:1000) Abcam Cat#ab11089; RRID:AB_2889189 Rat monoclonal anti-Ki67 (clone SolA15) antibody (mIHC; 1:1500) Thermo Fisher Scientific Cat#14-5698-82; RRID:AB_10854564 Mouse monoclonal anti-pan cytokeratin (PanCK) (clone AE1/AE3) antibody (mIHC; 1:200) Dako Cat#M3515; RRID:AB_2132885 Mouse monoclonal anti-alpha smooth muscle actin (aSMA) (clone 1A4) antibody (mIHC; 1:2000) Sigma-Aldrich Cat#A5228; RRID:AB_262054 Donkey anti-Rat IgG (H + L) Highly Cross- Adsorbed Secondary Antibody, Alexa FluorTM 594 Thermo Fisher Scientific Cat#A-21209; RRID:AB_2535795 Donkey anti-Rat IgG (H + L) Highly Cross- Adsorbed Secondary Antibody, Alexa FluorTM 488 Thermo Fisher Scientific Cat#A-21208; RRID:AB_2535794 Donkey anti-Rabbit IgG (H + L) Highly Cross- Adsorbed Secondary Antibody, Alexa FluorTM 568 Thermo Fisher Scientific Cat#A-10042; RRID:AB_2534017 Donkey anti-Rabbit IgG (H + L) Highly Cross- Adsorbed Secondary Antibody, Alexa FluorTM 488 Thermo Fisher Scientific Cat#A-21206; RRID:AB_2535792 Goat anti-Rat IgG (H + L) Cross-Adsorbed Secondary Antibody, Alexa FluorTM 568 Thermo Fisher Scientific Cat#A-11077; RRID:AB_2534121 Chemicals, peptides, and recombinant proteins Collagenase Sigma-Aldrich Cat#C7657 Paraformaldehyde Thermo Fisher Scientific Cat#28908 Penicillin streptomycin Gibco Cat#15140-122 (Continued on next page) 16 Cell Reports 43, 114837, October 22, 2024 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Ascorbic acid Sigma-Aldrich Cat#A4544 L-glutathione Alfa Aesar Cat#J62166 FBS Sigma-Aldrich Cat#F7524-500ML DMSO Chem Cruz Cat#358801 Triethanolamine Sigma-Aldrich Cat#90279 MetaPhor Agarose Lonza Cat#50181 Critical commercial assays Opal 6-plex detection kit Akoya Biosciences Cat#NEL871001KT Deposited data Deconvoluted raw light sheet imaging data This paper https://doi.org/10.6019/S-BIAD1136 Software and algorithms Fiji Schindelin et al.75 https://fiji.sc/ WEKA Segmentation plugin Arganda-Carreras et al.76 https://imagej.net/plugins/tws/ TrakEM2 Schindelin et al.75 https://imagej.net/plugins/trakem2/ Imaris FilamentTracer Bitplane ltd. RRID: SCR_007366 Tree Surveyor Short et al.77 Available from developers by request Qupath Bankhead et al.78 https://qupath.github.io/ MATLAB R2019a vol3D v2 script Woodford et al.79 https://se.mathworks.com/matlabcentral/ fileexchange/22940-vol3d-v2 3D branching simulation v. 3.21 This paper https://doi.org/10.17632/y8d52ygdxv.1 GraphPad Prism software package (version 10.1.2.) GraphPad RRID: SCR_002798 Article ll OPEN ACCESSEXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS Human tissue samples At the Turku University Hospital, normal breast tissue was obtained from donors undergoing breast reduction mammoplasty. Normal breast tissue adjacent to the tumor was collected in mastectomy due to pre-invasive carcinoma. Tissues were voluntarily donated upon written informed consent at Turku University Hospital, where the collection of clinical samples has been approved by the ethics committee of the the wellbeing services county of Southwest Finland (23/1801/2018) and was permitted by Turku University Hospital (T107/2018-1). The donors filled a questionnaire providing details of their reproductive history and hormonal medication. At Mediclinic in Oud-Heverlee (BE), human breast tissues were anonymously retrieved from healthy donors who underwent volun- tary cosmetic reduction surgery in the absence of any medical indication. All donors were adequately informed by third party clini- cians about the use of the tissue and provided written informed consent for the use of the resection material and publication of the experimental results. Limited and fully anonymous donor information (age, parity status, familial cancer risk) was obtained without any link to the donor’s medical record. All tissues from 17 donors were surgically dissected to obtain pieces with glandular tissue. For samples with high apparent mammographic density, a mild pre-digestion in collagenase was utilized before fixation [1.5h with 300 U/ml collagenase (Sigma, C7657)]. METHOD DETAILS Tissue clearing FUnGI tissue clearing and immunolabelling of tissue pieces was done following previously described protocol.18 All incubation steps were performed on a roller mixer. Tissue pieces were fixed in 4% PFA in +4C for 24 h, and washed and stored in PBS +1:100 peni- cillin/streptomycin (Gibco, 15140-122) at +4C until labeling. Fixed tissue pieces were labeled with DAPI (1:1000-1:3000) in PBS for 2 days in +4C. The areas containing TDLU structures were visualized with fluorescence microscopy, and chosen areas trimmed into pieces of approximately 2–5 mm in diameter. Tissues were washed in wash buffer [WB; 50 mg/mL ascorbic acid (Sigma, A4544), 0.05 mg/mL L-glutathione (Alfa Aesar, J62166) in PBS +0.1% Tween] for 30 min, and in WB1 [50 mg/mL ascorbic acid, 0.05 mg/mL L-glutathione in PBS +0.2% Tween +0.2% Triton X-100 (Sigma, T8787) + 0.2% SDS (Amresco, 0227) + 0.2% BSA (Biowest, P6154)] for 2–3 h in +4C, and incubated in anti-keratin 8 primary antibody [1:100 (Hybridoma Bank, TROMA-1)] diluted in WB2 (50 mg/mL ascorbic acid, 0.05 mg/mL L-glutathione in PBS +0.1% Triton X-100 + 0.02% SDS +0.2% BSA) for 2 days in +4C. After,Cell Reports 43, 114837, October 22, 2024 17 Article ll OPEN ACCESSpieces were washed 33 1 h in WB2, and incubated in anti-rat secondary antibody [1:400 (Thermo Fisher Scientific, Alexa Fluor 594)] diluted in WB2 for 36 h. Pieces were then washed again 3 3 1 h in WB2, after which they were moved to 25–30 mL of freshly pre- pared FUnGI solution [FUnGI base: 50% glycerol (Honeywell, 15675690), 2.5 M fructose (ACROS Organics, 1064642), 2.5 M urea (ACROS organics, 424585000), 10.6 mM Tris Base, 1 mM EDTA. Added prior to clearing: 50 mg/mL ascorbic acid, 0.05 mg/mL L-glutathione] and incubated for 2 h in room temperature. Finally, pieces weremoved to25–30mL of fresh FunGI solution and incu- bated at minimum o/n in +4C, or until mounting. iFLASH tissue clearing and immunolabelling were adapted from previously described protocol.21 All incubation steps were performed on a roller mixer. Tissue pieces were fixed with either 4% PFA in SDS permeabilization buffer (10% SDS in H2O, pH 7.4) o/n in +4 C or with 2% PFA in SDS permeabilization buffer o/n in room temperature. The next day, pieces were washed 3 3 5 min with PBS and moved to SDS permeabilization buffer for 1–2 days in +37C. Fixed tissue pieces were then transferred to 25-50mL of CUBIC1 [1.6 M urea, 5% Quadrol (Sigma, 122262), 15% Triton X-100, 25 mM NaCL (Fisher Scientific, S/3160/60) in ddH2O]80 and incubated in +37C for a minimum of 3 weeks, changing the buffer every 2–3 days. Next, the pieces were labeled with DAPI [1:1000-1:3000, (Life Technologies, D1306)] in PBS o/n at room temperature. The areas containing TDLU structures were visualized with fluorescence microscopy, and chosen areas cut into pieces of approximately 2–5 mm in diameter. Tissues were then permeabilized again using SDS permeabilization buffer for 2 days in +37C, changing the buffer between the days. Then, the tissues were washed 3 3 1 h in PBT (PBS +0.2% Triton X-100, pH 7.4) at room temperature and blocked with iFLASH blocking buffer [10% FBS (Sigma, 122262), 5% DMSO (Chem Cruz, 358801), 0.1% NaN3, 1% BSA in PBT] for 1 h at room temperature, and incubated in primary antibodies diluted in iFLASH blocking buffer for 3 days [anti-keratin 8 primary antibody 1:100 (Hybridoma Bank, TROMA-1); anti-keratin 14 primary antibody 1:100 (Bio- legend, 905301), anti-TAGLN antibody 1:50 (Abcam, ab155272)] in +37C. Tissues were washed 3 3 1h with PBT at room tem- perature, and pieces incubated in secondary antibodies [anti-rat secondary antibody 1:400 (Thermo Scientific, Alexa Fluor 488), anti-rabbit 568 secondary antibody 1:400 (Thermo Scientific, Alexa Fluor 568)] in +37C. Finally, the pieces were washed 3 3 1 h with PBT at room temperature, and transferred to 25-50mL of CUBIC2 [1.2 M sucrose (Millipore, 107651), 3.6M urea, 9% trie- thanolamine (Sigma, 90279), 0.1% Triton X-100 in ddH2O] for a minimum of 2 days in +37C, and moved to +4C until mounting and imaging. Labeled and cleared tissue pieces were imaged again to trim and position them optimally for the light sheet imaging excitation and emission light paths. Pieces were mounted in custom sample holders in 2% agarose (Lonza, 50181) diluted in milliQ water (FUnGI protocol) or CUBIC2 (iFLASH protocol), and let set at +4C. Finally, the sample holders were transferred to 25–30 mL fresh FUnGI or CUBIC2 respectively and incubated at minimum o/n in +4C, or until imaging. Histology and immunohistochemistry Tissue samples were fixed overnight at +4C in periodate–lysine–paraformaldehyde (PLP) buffer [1% paraformaldehyde (Thermo Scientific, 28908), 0.01M sodium periodate (Sigma, 311448), 0.075 M L-lysine (Alfa Aesar, 15404639) and 0.0375M P-buffer (0.081 M Na2HPO4 and 0.019M NaH2PO4; pH 7.4)]. After washing twice with P-buffer, samples were incubated in 30% sucrose in P-buffer for a minimum of two days. Within 14 days, samples were mounted in Tissue-Tek O.C.T. Compound (Sakura, 4583) on dry ice and cut into 8 mm sections. Frozen sections were H&E-labelled with standard procedures. For immunofluorescence labeling, the frozen sections were first thawed for 1 h at room temperature. The sections were blocked and permeabilized in permeabilization buffer (2% BSA, 0.1% Triton X-100 in PBS) for 30min at room temperature. Primary antibodies were incubated in 2%BSA in PBS overnight at +4C [anti-keratin 8 primary antibody 1:1000 (Hybridoma Bank, TROMA-1); anti-keratin 14 primary antibody 1:1000 (Biolegend, 905301)]. Then, the sec- tions were washed 3 3 10 min with PBS and incubated with the secondary antibodies, diluted in 2% BSA in PBS, for 1 h at room temperature [anti-rabbit secondary antibody 1:400 (Thermo Scientific, Alexa Fluor 488), anti-rat 568 secondary antibody 1:400 (Thermo Scientific, Alexa Fluor 568)]. Sections were then washed 3 3 10 min with PBS (1:1000 DAPI in the second wash) and then 5minwithmilliQ water. Finally, the sectionsweremounted under a #1.5 glass coverslip withMowiol (Calbiochem) supplemented with 2.5% 1,4-diazabicyclo[2.2.2]octane (DABCO, Sigma-Aldrich). Light sheet microscopy and deconvolution Samples were imaged with MSquared Aurora Airy Beam Light Sheet microscope using Aurora acquisition software version 0.5. The objective used was a Special Optics dipping objective which has a dipping medium refractive index (RI) dependent magnification 15.33-17.93 (NA 0.37–0.43) in immersion medium RI range 1.33–1.56. Alexa 488 secondary antibody was imaged with a 488 nm laser and emission filter at 500–540 nm and Alexa 594/Alexa 568 secondary antibodies with a 568 nm laser with long pass emission filter at 570 nm. 3D images were acquired by taking z-stacks with 400 nm spacing with a pixel size of 387 nm. Images were decon- voluted using Aurora Deconvolution software version 0.5 using point spread functions of fluorescent beads obtained from agarose embedded samples. Beads (Fluoro-max Fluorecent Polymer Microspheres, Thermo Scientific) were diluted 1:10 000 - 1:20 000 in 2% agarose (as described previously) and vortexed well for uniform distribution. Mixture was let settle for 10 min at room temperature to allow bub- bles to rise to the surface, and then mounted into an Eppendorf lid. The solidified pieces were then cut out. Additionally for FunGI protocol, bead samples were incubated in FUnGI o/n before imaging.18 Cell Reports 43, 114837, October 22, 2024 Article ll OPEN ACCESSConfocal imaging Cleared tissue samples were placed in 3.5mm glass-bottom dishes (Cellvis, D35-14-1-N) for imaging, additionally under a #1.5 coverslip covered in Vectashield (Bionordika, H-1000-10) if necessary to reach maximum imaging depth. Samples were imaged with the 3i Marianas CSU-W1 spinning disk (50mm pinholes) confocal microscope using SlideBook 6 software. The objective used was a 103 Zeiss Plan-Apochromat objective (NA 0.45). Signal was detected with the following solid-state lasers and emission filters; 405nm laser, 445/45nm; 488nm laser, 525/30nm; 561nm laser, 617/73nm. 3D images were acquired by taking z-stacks with a 2mm step size using a Photometrics Prime BSI sCMOS camera with a pixel size of 6.5 mm. Acquired images were 16bit with pixel dimensions 2048 3 2048. Tissue sections labeled with immunofluorescence were imaged at the same setup using the 20x Zeiss Plan-Apochromat (NA 0.8) and 63x Zeiss Plan-Apochromat (NA 1.4) objectives. Images were acquired by taking a z stack with 6 planes with the step size of 1 mm. Image segmentation and branching analysis Only images reliably showing a substantial fraction of the TDLU volume were chosen for branching analysis (n = 67). For morpho- metric measurements, segmentation of the 3D image was done using the Trainable WEKA Segmentation plugin.76 Images were pre- processed with background subtraction by rolling ball algorithm. A random forest pixel classifier was trained individually to each 3D image. 3–5 of the images in each image stack were used as training data, and a 200 tree random forest classifier was trained using Gaussian blur, mean,maximum,minimum,median, Laplacian, entropy, Sobel, difference of Gaussians, bilateral and Kuwahara filters as features. All segmentations were manually checked and fixed for mistakes. Quantitative analysis of the segmented TDLU structures was done by FilamentTracer function of Imaris 8.1.2 and 10.0.0 (Bitplane ltd.). After automatic placement of seed points using a minimum size of 25 mm estimated from the minimum branch diameter in the data, seed points weremanually checked and editedwhen necessary. Model surfacewasmanually thresholded tomatch the surface of the TDLU branches in the data. Branch diameters were calculated as the shortest distance from the distance map. Resulting fila- ment models were extensively manually checked for errors and corrected before quantitative measurements were exported. For morphometric measurements, FUnGI-cleared TDLUs showing the whole TDLU volume, and in which the connectedness of the branches could be reliably traced were selected for analyses concerning topological properties of the whole TDLU (n = 8). Filament models of TDLU structures were exported as TIFF stacks from Imaris and further analyzed using Tree Surveyor version 1.0.12.177 to extract parameters of subtree ancestry, network connectedness and the branching rotation angles. Rotation angles were extracted as the sharp (0–90⁰) angle between two planes defined by branches in two successive branching events. Localization of KRT14hi cells was observed from analyzed filament models in Imaris, and matched to corresponding subtree ancestry and network connectedness in Tree Surveyor. Total branch count of each TDLU was extracted from finalised filament models. TDLU volume was estimated by perpendicularly measuring the length and width of each structure and assuming a regular shape applying the formula for an ellipsoid’s volume; VTDLU = p 6 $length$width2Paraffin embedding of cleared tissue samples After light sheet imaging, cleared tissue pieces were stored in CUBIC2 in +4C for up to 4 months. For re-embedding into paraffin, pieces were manually extracted from the agarose mounting and remounted in HistoGel (Epredia, HG-4000-012) in desired orienta- tion. HistoGel-mounted tissue pieceswere re-embedded into paraffin blocks, fromwhich a 3mmcore containing the tissue of interest was taken,moved to a newmold,melted down and re-embedded into paraffin in a TMA formation. The constructed TMAblockswere then serial sectioned with a 4 mm section thickness up to 1.5mm of depth, or until glandular tissue was finished. Every 10th section was then H&E stained for 3D reconstruction and imaged using the Pannoramic P1000 slide scanner, and the remaining sections used for Opal 6-plex immunolabelling at planes of interest chosen from H&E. 3D reconstruction The 3D reconstruction for the TMA spot selected for visualization was done following a pipeline modified from.81 First, the stack was non-linearly aligned using the TrakEM2 plugin in Fiji/ImageJ v2.14.0.75 Following registration, acini were manually segmented in Qupath v0.5.1 using the brush tool.78 The acquired masks were then employed to discard background regions from the images in Python. Finally, the 3D reconstruction was created in MATLAB R2019a using the Vol3D v2 script.79 Multiplex staining Multiplex immunohistochemistry (mIHC) staining was performed using the BOND RX automated stainer (Leica Microsystems Ltd). For the multiplex mIHC the Opal 6-plex detection kit was used (Akoya Biosciences, NEL871001KT) and staining was conducted according to the manufacturer’s protocol. All reagents were prepared in advance and slides were loaded on the instrument. To begin the staining process, PFA-fixed and paraffin embedded tissue sections (4 mm) were incubated for 1h at 60C with Bond Dewax solution (Leica Bond, AR9222) to deparaffinize tissue sections. Antigen retrieval was achieved by immersing the sectionsCell Reports 43, 114837, October 22, 2024 19 Article ll OPEN ACCESSwith BOND Epitope Retrieval Solution 1 (Leica Bond, AR9961) for 20 min at 95C. Subsequent blocking of non-specific binding sites was followed for 15 min at RT. According to the primary antibody species, a different blocking buffer was used (for mouse and rabbit primary antibodies: Akoya Biosciences, ARD1001EA; for rat primary antibodies: Akoya Bioscience, 30024). Next, slides were incubated for 30 min at RT with primary antibodies [anti-p53 antibody 1:500 (Santa Cruz, sc-25268); anti-ER antibody 1:50 (Cell Signaling, 13258S); anti-KRT19 antibody 1:1000 (Thermo Fisher Scientific, MA5-12663); anti-FSP1/S100A4 antibody (Sigma, SAB5701199); anti-E-cadherin antibody 1:3000 (Abcam, ab76319); anti-KRT8 antibody 1:1000 (Sigma, MABT329); anti-CD3 anti- body 1:1000 (Abcam, ab11089); anti-Ki67 antibody 1:1500 (Thermo Fisher Scientific, 14-5698-82); anti-pan cytokeratin antibody 1:200 (Dako, M3515); anti-alpha smooth muscle actin antibody 1:2000 (Sigma, A5228)] which was diluted in TNB buffer (Akoya Biosciences, FP1020), followed by incubation with secondary antibodies ImmPRESS-HRP anti mouse/rabbit (Akoya Biosciences, ARH1001EA) or ImmPRESS- anti rat HRP (Akoya Biosciences, 30032) for 10 min at RT. Following this step, an Opal fluorophore working solution (Akoya Biosciences, NEL871001KT) was applied for 10 min at RT. All Opal reagents were diluted 1:150 in an amplification buffer (Akoya Biosciences, FP1609) except Opal 780, which was diluted to 1:50 in TNB buffer. To finish the cycle, primary and secondary antibodies were stripped with BOND Epitope Retrieval Solution 1 for 20 min at 95C. Subsequent staining rounds followed the same steps. Finally, nuclei were counterstained with DAPI (Akoya Biosciences, FP1490) for 5 min at RT before mounting the stained tissue with mounting medium (Invitrogen, P36961). Stained sections were imaged using a slide scanner (Akoya Vector Polaris). Multiplex image analysis The branch counts of singular TDLUs were determined from consecutive H&E sections by adding up all distinguishable acini per structure. Matching TDLUs were then located in multiplex imaging and cropped out for individual analysis in Fiji. Then, relevant epithelial area was segmented; ER to luminal epithelium (KRT19), p63 to basal epithelium (aSMA), Ki67 to epithelium (pan-keratin), CD3 to within [epithelium] or outside of luminal epithelium [intra-TDLU stroma] (KRT8) and S100A4 to within [epithelium] or outside of epithelium [intra-TDLU stroma] (CDH1). The resultingmaskwas then applied to signal of interest, fromwhich either its area (p63, CD3, S100A4) or number of nuclei (ER, Ki67) were determined. For stromal signal analysis, intra-TDLU was determined by connecting pe- ripheral acini, and then removing the signal of interest from within the epithelium. Mathematical modeling BARW simulations were implemented as in42 including branch self-repulsion.45 The simulated networks will grow indefinitely unless growth is somehow limited. We tested two mechanisms for network growth limiting: 1) Growth was limited to a predefined volume, and 2) a model in which branches have branching potential that wanes in older generations and eventually stops the branching. 1000 simulations were run for both models to generate the dataset. (1) Volume-limited model At each simulation step each active branch can do one of three possibilities: (1) it can grow by one more step, (2) it can bifurcate to two child branches or (3) it can terminate (i.e., stop growing and become inactive) if it meets an existing branch or limits of the simu- lation volume. Each branch has a propagation direction unit vector bkðtÞ pointing to the direction of branch growth at time t. At each simulation step a random number in [0,1] is drawn and bifurcation happens if the number is smaller than the bifurcation probability. Otherwise the branch grows forward one step in direction bk. Similar to,42 at each step the branching angle is diffused by a random angle in range [0, 15⁰] (angle range determined as in46) and a rotation angle in range [0, 360⁰] to prevent unnatural line-straight branches. If bifurcation happens, two bifurcation angles b are drawn from a Gaussian distribution fitted to experimental bifurcation angle distribution (Figure S5H), and two angles q are drawn from a uniform angle distribution of [0, 360⁰] to represent the rotational directions of the new branches around the parent branch (Figure S5I). An active branch terminates if its tip gets within termination distance rT from an existing branch network point. Simulation starts with one active tip starting to grow with bk = bex. Simulation runs until there are no active tips left. 3D vector rotations are done by the Rodrigues’ vector rotation.82 vrot = v cos a + ðba 3 vÞsin a+ baðba$vÞð1  cos aÞ (Equation 1) where v is the original vector to be rotated by angle a around the axis defined by the unit vector ba, and vrot is the final rotated vector. Direction vectors bk of the new branches in 3D space are calculated using the following algorithm: first a reference unit vector defining a direction perpendicular to parent branch direction is calculated as bu = bk 3 bey , where bey is the standard unit vector along y axis. Vector bu is then rotated by the rotation angle q around parent branch’s bk vector using Equation 1 to get a randomly in [0.360⁰] oriented rotation axis ba for the final rotation with the branching angle b. The direction of the child branch is then calculated by rotating the parent branch’s direction vector bk by branching angle b around the axis ba again using Equation 1. Simulation parameter values were derived from the experimental data as follows and they were scaled to represent units of mm. Branch diameter value of 40 mmwas used which was calculated by adding an extra 5 mm for the basal cell layer to the experimentally observed 35 mm branch diameter of the luminal cell layer. Bifurcation probability was calculated to be 0.15 by dividing the simulation step size 10 mm with the experimentally observed average branch length of 71 mm.20 Cell Reports 43, 114837, October 22, 2024 Article ll OPEN ACCESSBranch repulsion causes an additional displacement to the direction vector bk. Every active tip experiences repulsion from all other existing network points that are within repulsion radius in the front half-space of the tip. Displacement due to repulsion d ! R was calcu- lated as d ! R =  fR bR , where fR is the repulsion strength and bR unit vector pointing from the examined tip to the mass center of all network points within repulsion radius in the front half-space. d ! R was added to bk and the result was normalized to give direction unit vector for the next step. fR was optimized by finding the best match between the branch count distributions of the simulated and experimentally observed networks. Termination distance was determined from the experimental data by finding terminated branches that pointed toward another branch and measuring the gap from the tip of such branch to the surface of the second branch. The average gap in such cases was 18 mm (SD 11 mm) (Figure S5J). Since the calculations in the model are done with branch paths that have thickness of zero (solid lines in Figure 5A), the termination distance was set equal to 38 mm consisting of the gap plus half the branch diameter (Figure 5A). To prevent immediate termination by the parent and sibling branch at the branching event, termination events only in the front half-space of the active tip were taken into account like in.46 As observed previously,46 simulation produces realistic networks only if immediate re-branching after a branching event is prevented and therefore a branching free zone of radius equal to branch diameter was estab- lished around each branching event to prevent immediate re-branching which sets a requirement that branches shorter than their width are not mature enough to branch again (Figure 5A).46 Self-repulsion from close by branches45 was used and the value of repul- sion strength was optimized to 0.2 by finding the best match between the branch count distributions of the simulated and experimen- tally observed networks (Figure S7). Repulsion distance was determined from the experimental data by finding branches that had been growing toward another branch but making a bend when getting close. The distance from the bend to the surface of the second branch was on average 35 mm (SD 11 mm) (Figure S5J). The repulsion distance was set equal 55 mmconsisting of the gap plus half the branch diameter again. If a rigid predefined simulation volume box is used in the volume-limitedmodel to limit the growth, some networks that start to grow at an angle with respect to the box can terminate prematurely by colliding to the sides of the simulation volume. There we did not confine the simulation to a predetermined box, but instead the network was allowed to grow in the direction it chooses. Simulation volume was controlled by terminating any branch that wanders further in length from the origin (maximum length) or wider from the general direction of the network (maximum width) than allowed by the limits (Figure 5B). General direction of the network was deter- mined as a vector from origin to centroid of the set of current network points. Maximum length was set equal to experimentally observed average TDLU length of 460 mm and maximum width to 310 mm to create a limiting volume of typical TDLU size and having an aspect ratio of 1.5. In effect this creates a limiting volume that follows along the growing network. (2) Model of gradually waning branching potential Mechanics of the simulation were otherwise identical to the volume-limited model, but in the branching potential waning model we modified the BARWmechanism so that when a branch stops growing and goes to the bifurcation state, instead of always bifurcating it makes a decision to bifurcate or self-terminate. Probability to self-terminate increases with branching level so that the probability to branch follows equation pBðLÞ = qL, where L is the branching level and q is a waning rate constant. This describes a fractional decline of branching potential by increasing branching level. Setting a requirement that pBð0Þ = 1; we estimated the waning rate q = 0:867 from a fit to the experimentally observed distribution of KRT14hi cells (Figure 5C) in different levels. In the simulation, when a branch reaches the end of elongation state (determined by the bifurcation probability parameter as in the volume-limited model), then a random number R in [0,1] is drawn, and if R