UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modeling milk production in the lactation period and the effect of feeding frequency on milk production Faulkner, Katharine R. 2019

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2019_november_faulkner_katharine.pdf [ 901.04kB ]
JSON: 24-1.0380624.json
JSON-LD: 24-1.0380624-ld.json
RDF/XML (Pretty): 24-1.0380624-rdf.xml
RDF/JSON: 24-1.0380624-rdf.json
Turtle: 24-1.0380624-turtle.txt
N-Triples: 24-1.0380624-rdf-ntriples.txt
Original Record: 24-1.0380624-source.json
Full Text

Full Text

Modeling milk production in the lactation periodand the effect of feeding frequency on milkproductionbyKatharine R. FaulknerA THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THEREQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2019c© Katharine R. Faulkner, 2019The following individuals certify that they have read, and recommend to the Facultyof Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled:Modeling milk production in the lactation period and the effect of feeding frequencyon milk productionsubmitted by Katharine R. Faulkner in partial fulfillment of the requirements for thedegree of Master of Science in Mathematics.Examining Committee:Eric Cytrynbaum, MathematicsSupervisorLeah Edelstein-Keshet, MathematicsSupervisory Committee MemberNeil Balmforth, MathematicsSupervisory Committee MemberAdditional ExaminerAdditional Supervisory Committee Members:Supervisory Committee MemberSupervisory Committee MemberiiAbstractHuman milk production is controlled by a variety of internal and external factors,including hormones, neurons, suckling stimulus and milk removal. One method forincreasing milk production suggested to mothers who want to produce more milk iscluster feeding: splitting one feeding into a cluster of multiple feedings. Phenomeno-logical models have been proposed to describe average weekly milk production ratesin dairy cattle, but these models do not take into account the effect of the milkremoval schedule used. In this thesis, two ordinary differential equation models fordescribing milk production and milk removal are presented: Model 1, which assumesa linear rate of milk removal during feeding, and Model 2, which uses physical modelsof fluid flow through a system of ducts to estimate the rate of milk removal. Thesemodels are qualitatively very similar, but Model 2 allows for examination of differ-ences in milk dynamics between alveoli at different depths in the mammary gland.24 hour milk production was then predicted for each model version using variouscluster feeding schedules. Feeding schedules which had the most frequent and regu-lar feedings elicited small increases in milk production during a day compared to lessfrequent or regular schedules. This small increase in milk production suggests thatcluster feeding may not be an effective method of increasing production, as other fac-tors that decrease production, like stress, may override this small effect. More workneeds to be done to test these models in order to determine an effective method forincreasing milk production.iiiLay SummaryMathematical models of physiological systems allow researchers to gain insightsinto the system being studied through mathematical analysis, as well as performexperiments via theoretical simulations. In this thesis, two models for human milkproduction are presented and tested. One method for increasing milk productionsuggested to mothers who want to produce more milk is cluster feeding, splitting onefeeding into a cluster of multiple feedings, so these models were tested using variouscluster feeding schedules. Feeding schedules which had more frequent feedings elicitedsmall increases in milk production compared to less frequent schedules. This smallincrease in milk production suggests that cluster feeding may not be an effectivemethod of increasing production, as other factors that decrease production, likestress, may override this small effect. More work needs to be done to test thesemodels in order to determine an effective method for increasing milk production.ivPrefaceThis thesis is original, unpublished work by the author, K.R. Faulkner.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 Model 1: Constant milk removal . . . . . . . . . . . . . . . . . . . . . 42.1.1 No milk age signal . . . . . . . . . . . . . . . . . . . . . . . . 82.1.2 Milk age affecting cell quiescence . . . . . . . . . . . . . . . . 92.1.3 Milk age affecting alveolus size . . . . . . . . . . . . . . . . . . 92.1.4 Milk age affecting cell quiescence and alveolus size . . . . . . . 102.2 Model 2: Pressure-based milk removal . . . . . . . . . . . . . . . . . 102.3 Parameter Choices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.1 Model 1: Constant milk removal . . . . . . . . . . . . . . . . . . . . . 163.1.1 No milk age signal . . . . . . . . . . . . . . . . . . . . . . . . 163.1.2 Milk age affecting cell quiescence . . . . . . . . . . . . . . . . 193.1.3 Milk age affecting alveolus size . . . . . . . . . . . . . . . . . . 213.1.4 Milk age affecting cell quiescence and alveolus size . . . . . . . 233.1.5 Cluster feeding with constant milk removal . . . . . . . . . . . 253.2 Model 2: Pressure-based milk removal . . . . . . . . . . . . . . . . . 313.2.1 Flow out of a single alveolus . . . . . . . . . . . . . . . . . . . 313.2.2 Flow out of multiple alveoli . . . . . . . . . . . . . . . . . . . 354 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.1 Model 1: Constant milk removal . . . . . . . . . . . . . . . . . . . . . 414.2 Model 2: Pressure-based milk removal . . . . . . . . . . . . . . . . . 434.3 Increasing milk production . . . . . . . . . . . . . . . . . . . . . . . . 444.4 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46A Parameter Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51viList of TablesA.1 Milk production model physical constants. . . . . . . . . . . . . . . . 51A.2 Milk production model parameters. . . . . . . . . . . . . . . . . . . . 51viiList of Figures2.1 Diagram of key variables in milk production with locations in body. . 52.2 Diagram of pressures and duct parameters in branching alveoli network. 133.1 Simulated milk for 24 hours of feedings with Model 1 and no milk aging. 163.2 Phase lines for Model 1, constant milk removal with no milk aging,and the two system states: (a) between feedings (no milk removal)and (b) during feeding (sustained removal). . . . . . . . . . . . . . . . 183.3 Simulated milk curve for the end of feeding, as predicted by Model 1,constant milk removal without milk aging. . . . . . . . . . . . . . . . 193.4 Simulated milk curve for 24 hours of feeding one hour in every four, aspredicted by Model 1, constant milk removal with aged milk affectingcell quiescence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.5 Simulated milk curve for the end of feeding, as predicted by Model 1,constant milk removal with aged milk affecting cell quiescence. . . . . 213.6 Simulated milk curve for 24 hours of feeding one hour in every four,with Model 1, constant milk removal with aged milk affecting alveolussize. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.7 Simulated milk curve for the end of feeding, as predicted by Model 1,constant milk removal with aged milk affecting alveolus size. . . . . . 233.8 Simulated milk curve for 24 hours of feeding one hour in every four,with Model 1, constant milk removal with aged milk affecting cellquiescence and alveolus size. . . . . . . . . . . . . . . . . . . . . . . . 243.9 Simulated milk curve for the end of feeding, as predicted by Model1, constant milk removal with aged milk affecting cell quiescence andalveolus size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.10 Total milk produced in 24 hours for each of the four versions of Model1, and six different feeding schedules. . . . . . . . . . . . . . . . . . . 263.11 Total milk produced in 24 hours for each of the four versions of Model1, with varied feeding schedules parameterized by the time betweenfeedings in a single cluster (tb) and the number of feedings in a cluster(Nc). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28viii3.12 Total milk produced in 24 hours as predicted by Model 1 with agedmilk affecting cell quiescence and alveolus size, and with varied feed-ing schedules parameterized by the time between feedings in a singlecluster (tb) and the number of feedings in a cluster (Nc). . . . . . . . 293.13 Simulated milk curve for 24 hours with Model 1, constant milk removalwith aged milk affecting cell quiescence and alveolus size. . . . . . . . 303.14 Simulated milk curve for 24 hours of feeding one hour in every four,with Model 2, pressure-based milk removal with a single alveolus. . . 323.15 Total milk produced in 24 hours for each of the four versions of Model2 in a single alveolus, with varied feeding schedules parameterized bythe time between feedings in a single cluster (tb) and the number offeedings in a cluster (Nc). . . . . . . . . . . . . . . . . . . . . . . . . 343.16 Simulated milk curve for 24 hours of feeding one hour in every four,with Model 2, pressure-based milk removal with three alveoli. . . . . 363.17 Total milk produced in 24 hours for each of the four versions of Model2 in three alveoli, with varied feeding schedules parameterized by thetime between feedings in a single cluster (tb) and the number of feed-ings in a cluster (Nc). . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.18 Simulated milk curve for 24 hours of feeding one hour in every four,with Model 2, pressure-based milk removal with three alveoli. . . . . 39ix1. IntroductionBreastmilk is considered the best source of nutrition for infants and lactationis the process by which milk is produced in and secreted from the mammary glandin the breast (Buckley 2015; Amis 2015). In developing countries, lack of maternalnutrition lead to problems with lactation, but low milk production has become anissue for many mothers in developed countries despite adequate nutrition (Walker1997; Kent et al. 2012; Hurst 2007). Low milk production is especially prevalentin mothers of preterm infants, who must initiate milk production and continue toproduce milk throughout the stress from complications of preterm birth (Fewtrellet al. 2016). In these cases, mothers have access to enough food and supplements tobe able to produce milk at full capacity, but struggle to do so. There must then befactors dictating the amount of milk produced other than simply nutrition. Mothersand doctors in the developed world are then interested in methods of increasing milkproduction, and and such methods continue to be extensively studied (Daly et al.1996; Dewey 2001; Kent et al. 2012; Marasco 2015; Fewtrell et al. 2016).There are many methods that have been suggested for increasing milk produc-tion, such as herbal or pharmacological treatments to increase prolactin levels (thehormone which stimulates milk production), or increasing skin-to-skin contact be-tween the mother and infant (Kent et al. 2012). One method of increasing productionthat has been suggested is through manipulation of the feeding schedule—feedingmore frequently or clustering feedings may increase milk production (Kent et al.2012; Daly et al. 1996). This type of feeding schedule reduces the between breast-feeding, possibly allowing for more drainage of the breast, both of which are thoughtto improve milk production (Kent et al. 2012). To what extent these changes in feed-ing schedule might help increase production, and why increased drainage and fewerextended periods between feedings increase milk production is not fully understood.The goal of this thesis is then to build a model of lactation that explores possiblemechanisms behind increased milk production due to changes in the feeding schedule,and to determine how the feeding schedule can be altered to increase production.1Although milk production has been modeled before, it has been largely phe-nomenological models of dairy cow lactation that describe how milk productionchanges over the course of the entire lactation period, from initiation to weaning(Wood 1968; Dijkstra et al. 2010; Ferreira et al. 2015; Dag et al. 2005). These studieshave included fitting milk production data to very simple differential equation orstatistical models, and explored effects like season of lactation initiation on overallmilk production. These models also are made to explain changes in milk productionover the course of weeks or months, and cannot explain changes that occur on smallertime scales. To explore different types of feeding schedules, a model must be able toreflect changes in milk production that happen over hours or minutes.There is also a model of the fluid dynamics of milk removal in a human breast(Negin Mortazavi et al. 2016). In this model, the entire mammary gland was sim-ulated, and suckling patterns within a single feeding studied. This model does notinclude how milk is produced or how that production is regulated. Lactation is acomplicated system that is regulated by multiple different types of factors, includingneurological signals, many hormones and external physical stimuli. In order to studymilk production, these methods of regulation must be included.Milk production is initiated and promoted by the hormone prolactin, produc-tion of which is controlled and inhibited by the neurotransmitter dopamine secretedfrom specific neurons in the brain (Grattan 2015; Egli et al. 2010). Regulation ofprolactin secretion thus goes through these neurons and any regulators that promoteprolactin must inhibit TIDA neuron release of dopamine in order to do so. Infantsuckling is an example of a neurological signal that promotes prolactin—as the infantsuckles on the nipple, the pressure initiates a neurological signal through activationof mechanosensors (Tay et al. 1996). The hormone oxytocin also plays a role in lac-tation, but not in production—oxytocin causes myoepithelial cells (a thin layer ofmuscle cells) surrounding the alveoli to contract, pushing milk out of the breast whensuckling begins (Tay et al. 1996). These signals are the basic components of lactationthat play an important role in milk production. These components, as well as otherfactors included in the presented models will be described in more detail in the nextsection.2Mechanisms that will be explored are based on possible effects of milk ag-ing (Daly et al. 1993). Because prolonged periods between feedings and incompletedrainage of the breast inhibit milk production, the mammary gland may have someway of measuring and responding to milk aging as it accumulates in the breast. It hasbeen proposed that there is a signal that accumulates in milk as it ages and that themilk-producing cells in the mammary gland may respond to this signal (Daly et al.1993; Schmitt-Ney et al. 2014; Cox et al. 1996; Wilde and Peaker 1990). One culpritfor this signal is leukemia inhibiting factor (LIF), which has been implicated in cellcycle signalling pathways in the mammary gland of mice (Schere-Levy et al. 2003;Watson 2006). Whatever the signal may be in humans, it is proposed to have anautocrine (local) signalling effects rather than neurological or hormonal effects (Dalyet al. 1996; Kent et al. 2012; Wilde and Peaker 1990). Two possible local effects areconsidered: milk aging causing mammary cells to become quiescent (not active andwith the cell cycle paused), and aging causing apoptosis (programmed cell death) ofthe mammary cells.In this thesis, two main models are presented. Model 1 explores the physiologi-cal mechanisms behind milk production, specifically the local effects of milk aging onproduction. Four versions of Model 1 are compared in order to analyze different pos-sible effects that milk aging may have on mammary gland cells. This model considersa single alveolus in the mammary gland that is producing and storing milk. Feedingis modeled with a prescribed, constant rate of milk removal during the tested feedingschedule. Model 2 takes the final version of Model 1 and adds on a fluid dynamicsbased model of milk removal that takes into account some of the geometry of themammary gland. With this model, more alveoli can be simulated and the differencesbetween those alveoli explored.32. ModelsTwo main categories of models were explored: (1) constant flow rate of milkduring feeding, and (2) milk flow rate, derived from Poiseuille’s Equation for fluidflow at low Reynold’s number and pressure balance in the alveoli. Both of thesecases used the same general model to describe the relationship of the hormonal andneurological effects on milk production, as described below in Model 1. Model 2 usesthe same foundation model, and uses principles from fluid physics to allow for theconsideration of more complex geometries.2.1. Model 1: Constant milk removalMilk production is controlled by multiple factors, at various time scales. Milkis produced over the course of hours, and the production rate changes over thecourse of days (Dijkstra et al. 2010). On the other hand, the tuberoinfundibulardopamine (TIDA) neurons, which are the main neurological source of control formilk production, respond within milliseconds to changes in stimuli. The fast timescale variables include blood hormone levels (prolactin and oxytocin), which respondin minutes, and neural stimuli (dopamine, the suckling signal, and the stretch signal),which respond in milliseconds (Cox et al. 1996). The slow time scale variables includemilk production, the fraction of alveolar cells that are active and the amount ofalveolar cells, which respond in hours. Figure 2.1 shows the locations of release andaction for some of the variables involved in milk production.4Figure 2.1: Diagram of key variables in milk production with locations in body.Prolactin is produced in the anterior pituitary gland of the brain, and is regulatedlocally by dopamine released from the TIDA neurons in the pituitary gland. Prolactinis then released into the bloodstream to act at the mammary gland. Stimuli fromthe breast, including suckling and stretch of the alveoli, travel from the breast to thebrain via neurons and the spinal cord.Fast time scale variablesIn this model P represents the level of prolactin, the hormone directly controllingmilk production, in the bloodstream. This variable is unitless, corresponding to pro-lactin concentration relative to the maximum prolactin concentration. Previous mod-els have described prolactin dynamics in non-lactating humans and rats (Rattanakuland Lenbury 2009; Bertram et al. 2006). While some models account for regulationof prolactin by the hormone thyrotropin-releasing hormone (TRH), others considerTRH merely an intermediary step between dopamine regulation and prolactin secre-tion that can be ignored (Bertram et al. 2006). In Model 1, factors that could onlyaffect TRH were not considered, so for simplicity TRH was not included as a sepa-rate variable. Instead the parameters in Equation 1.1a are assumed to be functionsof TRH, and all regulation of prolactin secretion is through dopamine.Prolactin secretion from the posterior pituitary is under perpetual inhibition5from the TIDA neurons in the pituitary, which release dopamine according to theneurons’ firing rate, D (normalized to its maximum firing rate so D is unitless).Prolactin is then produced at a rate that decreases as the TIDA neuron firing rateincreases. Prolactin also is reabsorbed from the bloodstream at a rate proportional toitself. Similar to the model made by Bertram et al in 2006, the differential equationfor prolactin level isdPdt=kP1kP2 +D− kP3P. (1.1a)The constants kP1 and kP2 describe dopamine-regulated prolactin production so thatkP1/kP2 is the maximum prolactin production rate and kP2 is the relative TIDA firingrate at which prolactin is produced at half of its maximum rate. The rate of prolactinremoval is then kP3. In this model, it is assumed that prolactin levels in the bloodequilibriate quickly relative to milk production and feeding (Tay et al. 1996). Theright hand side of Equation 1.1a is set to zero and prolactin is at a quasi-steady statedescribed byP =kP1kP3(kP2 +D). (1.1b)Milk removal is also known to influence milk production, and accumulation of milkcould influence production either through stretch on the alveolus cortex or accumu-lation of a milk aging factor (Wilde and Peaker 1990). The effect of milk aging asa local signal will be discussed in a later section, but the stretching of the alveoluswas considered as a potential neurological signal that affects TIDA neuron firing inthis model. Additionally, the mechanical suckling on the nipple is known to influ-ence prolactin production, and would also do so through the TIDA neurons (Tayet al. 1996). These signals—stretch and suckling—would be sent via neurons fromthe breast to the TIDA neurons in the brain, so S is the relative firing rate of neuronscarrying the stretch signal and U is the relative firing rate of neurons carrying thesuckling signal (each normalized to their respective maximum firing rates, so S andU are unitless). Stretch of the alveoli indirectly decreases prolactin production andtherefore promotes TIDA neuron firing. Suckling on the nipple indirectly increasesprolactin production, through inhibition of TIDA neuron firing. These two stimuli,stretch and suckling, are assumed to affect TIDA neuron firing independently, with6each sufficient to affect TIDA firing on their own. The differential equation for D,relative TIDA neuron firing, isdDdt= kD1(1kD2 + U+SkD3 + S)+ kD4 − kD5D. (1.2a)Action potentials in the TIDA neurons are then produced at a rate inversely pro-portional to U and at an additional rate that saturates with S. Because prolactinproduction is under perpetual inhibition from the TIDA neurons, these neurons musthave a constant background action potential production rate kD4. In this model, itis assumed that the firing rate in the TIDA neurons equilibriates quickly relative tomilk production and feeding, so the right hand side of 1.2a is set to zero and thisfiring rate is at a quasi-steady state described byD =kD1kD5(1kD2 + U+SkD3 + S)+kD4kD5. (1.2b)The firing rate for the suckling signal, U , is produced only when the infant is sucklingand is assumed to decay away rapidly after suckling stops. The equation for thesuckling signal firing rate is thenU = kUf(t), (1.3)where f(t) is a characteristic function for the feeding schedule, with f = 1 duringfeeding and f = 0 otherwise. The parameter kU is the fixed firing rate during suckling.The stretch signal, S, is a saturating function of stress on the membrane, σ,S =kS1σkS2 + σ. (1.4)For the constant milk removal case, the stress was assumed to be directly related tothe amount of milk in the alveolus, N , with σ = kσN .Finally, during feeding, high oxytocin levels in the blood, ω, causes the myoep-ithelial cells around the alveoli to contract, squeezing milk out of the mammary glandby generating a pressure on the alveoli whose role is described later (Tay et al. 1996).As with prolactin, ω is normalized to the maximum oxytocin concentration level andω is unitless. The suckling signal causes this rise in oxytocin, and the differentialequation for oxytocin isdωdt= kω1U − kω2ω. (1.5a)7The suckling signal firing rate produces oxytocin with a constant rate kω1 and decaysaway with constant rate kω2. In this model, it is assumed that oxytocin equilibriatesquickly relative to milk production and feeding, so the right hand side of Equation1.5a is set to zero and oxytocin is at a quasi-steady state described byω =kω1kω2U =kω1kω2kUf(t). (1.5b)Slow time scale variablesFour variations of Model 1 were considered in order to explore different possiblemechanisms for the effect of aged milk on production. Because milk aging has beenfound to be an autocrine (local) signal, it must affect the cells producing milk di-rectly (Wilde and Peaker 1990; Schmitt-Ney et al. 2014). Two possible age signalswhich could affect milk production are leukemia-inhibiting factor (LIF) and sero-tonin, both of which are local signals in the mammary gland which have been shownto be important for involution, the termination of milk production after feeding ends(Schere-Levy et al. 2003; Horseman and Collier 2013).Two possible mechanisms for such an age signal were considered: the age signaldeactivating alveolar cells, with the cells left occupying space in the alveolus and notproducing milk, but capable of being recruited to produce milk again; and the agesignal causing apoptosis of alveolar cells, which are then removed or shed, decreasingthe size of the alveolus. Each combination of these two mechanisms was explored,including the case where neither exists and milk aging has no effect.No milk age signalWith no age signal, all molecules in the milk were treated as the same andtracked collectively as N , the number of moles of “dry” (non-water) milk moleculesin the alveolus. This milk is produced at a rate that depends on both the serumprolactin level, P , and the number of productive cells surrounding the alveolus,which is related to the reference (non-stretched) radius of the alveolus, r0, and theproportion of alveolar cells that are active, C:dNdt= kNPCr20 −NQV− klN, (2.1a)8where V = 4pir30/3 is the total volume of the alveolus and Q is the volumetric flowrate out of the alveolus (both constant in this case, with Q/V = kFl). The milk, Nis removed with feeding at rate NQ/V . Milk also leaks out of the alveolus, betweengaps in alveolar cells, or is reabsorbed back into the surrounding cells with constantrate kl.Milk age affecting cell quiescenceIn the case where aged milk does have an effect, it must be tracked separatelyfrom new milk as G, the amount of aged milk in moles. New milk, N , is convertedto aged milk at some rate kG, thus modifying the equation for N to include thisconversion:dNdt= kNPCr20 −NQV− klN − kGN, (2.1b)Aged milk has the same rate of removal and leakage as new milk:dGdt= −GQV− klG+ kGN. (2.2)If the age signal promotes cell quiescence, but not cell removal, then the number ofalveolar cells (and therefore the reference radius) would remain constant. The variableaffected would then be the proportion of these cells which are actively producing milk,C. In the absence of an age signal, all cells would be recruited to produce milk, withprolactin acting as an activator of quiescent cells (Horseman and Collier 2013). Asage signal accumulates, it would then deactivate cells, decreasing C. The differentialequation for C is thendCdt= kC1P (1− C)− kC2CG, (2.3)where prolactin activates cells at rate kC1, and the age signal deactivates cells at ratekC2. In this case, the alveolus reference radius, r0 remains constant.Milk age affecting alveolus sizeIf the age signal promotes cell removal, then the number of cells surrounding thealveoli would change, but all the cells would be actively producing milk. The numberof cells surrounding the alveolus is proportional to its reference surface area, 4pir20.The number of cells, and therefore the surface area, is assumed to grow logisticallywith a rate constant that depends on prolactin, leading to the following differential9equation for the radius:dr0dt= (kr1 + kr2P )r0(r20,max − r20)− kr3r0G. (2.4)The cells have a background growth rate of kr1, with a “carrying capacity” or maxi-mum radius r0,max, and prolactin increases the growth rate by kr2P . The age signal,G, accelerates apoptosis and removal of cells, decreasing the radius with constant ratekr3. In this case, the proportion of alveolar cells that are active, C, would remainconstant.Milk age affecting cell quiescence and alveolus sizeIf both of the described mechanisms for age signal action are in effect, thenneither C nor r0 would be constant. To test these two mechanisms together, bothequations 2.3 and 2.4 were used.2.2. Model 2: Pressure-based milk removalIn order to explore the effects of including multiple alveoli in a branching struc-ture with this model, as well as to test the hypothesis that differences in milk removalbetween alveoli in the mammary gland affect milk production, the flow rate out ofeach alveolus must be determined. These flow rates are dependent on the geometryof the branching, therefore principles from fluid physics are needed to describe therelationship. These principles were applied first to a simple structure with a singlealveolus, then to larger branching structures with multiple alveoli.Flow out of a single alveolus.To determine the flow rate of milk from a single alveolus, through a duct, and outof the nipple, we invoke Poiseuille’s law. According to this law, at low Reynold’snumber, the flow rate, Q, of fluid through a single pipe isQ =pir48νl∆P, (3.1)where r is the radius of the pipe, l is the length of the pipe, ν is the viscosity ofthe fluid and ∆P is the change in pressure from one end of the pipe to the other.According to one study, viscosity of breast milk ranges between 1.2 × 10−3 and2.5× 10−3 Pa·s (Sunaric´ et al. 2016). Various duct lengths and radii were tested.10For a single alveolus, the change in pressure across a single pipe is the differencebetween the pressure from the alveolus, PA, and the pressure at the nipple from theinfant suckling, PB. The pressure from the infant is prescribed and various values ofthis pressure were tested. The pressure from the alveolus can be derived from forcebalance on the alveolus cortex, assuming that the volume of the alveolus equilibriatesquickly so that these forces are in balance. This force balance yields the followingpressure in the alveolus:PA =KhrA(r2Ar20− 1)+ Pω, rA ≥ r0Pω, rA < r0,(3.2)where rA is the stretched radius of the alveolus, r0 is again the non-stretched refer-ence radius, and Pω is some prescribed pressure exerted by the myoepithelial cellssurrounding the alveolus as they contract in response to oxytocin during feeding(Newton and Newton 1948). In order to find the stretched radius, rA, the quasi-steady-state volume equation for a single cell (Jiang and Sun 2013) was modified foran alveolus to include outflow of milk as follows:dVAdt= −β (∆P −∆Π)−Q = 0. (3.3)Here, ∆P is the difference in hydrostatic pressure between the interior and exteriorof the alveolus (Pin − Pex), ∆Π is the difference in osmotic pressure between theinterior and exterior of the alveolus (Πin−Πex), and β is a parameter describing thevolume of water moved across the membrane per unit pressure change. The differencein hydrostatic pressure, ∆P , is the difference between the alveolus interior pressure,PA (Equation 3.2), and the contraction pressure, Pω. The difference is then∆P = PA − Pω =KhrA(r2Ar20− 1), rA ≥ r00, rA < r0.(3.4)The osmotic pressure, Π, is related to the concentration of milk in the associatedregion. The osmotic pressure in the cells surrounding the alveolus, Πex, is assumedto be constant at pressure ΠC . The difference in osmotic pressures is then∆Π =NVART − ΠC = 3N4pir3ART − ΠC , (3.5)where N is the moles of milk inside the alveolus, VA is the volume of the alveolus, R11is the ideal gas constant, and T is body temperature. The flow out of the alveolus, Q,is related to the difference in pressure between the alveolus, PA, and the end of theduct leading away from the alveolus. In the case with a single alveolus, this pressureis the pressure from the infant suckling, PB. Using Poiseuille’s Law, Equation 3.1,the flow rate out of the alveolus isQ =pir4d8νld(PA − PB) =pir4d8νld(KhrA(r2Ar20− 1)+ Pω − PB), rA ≥ r0pir4d8νld(Pω − PB), rA < r0,(3.6)where rd and ld are the radius and length of the duct between the alveolus andthe nipple. We can now use equations 3.4-3.6 to write dVA/dt in Equation 3.3 as afunction of rA, as follows:0 =−βKhrA(r2Ar20− 1)+ β(3N4pir3ART − ΠC)− pir4d8νld(KhrA(r2Ar20− 1)+ Pω − PB), rA ≥ r0β(3N4pir3ART − ΠC)− pir4d8νld(Pω − PB), rA < r0.(3.7)This equation can then be solved at each time step for rA using a non-linear solver(fsolve in Matlab 2019a, Mathworks). The flow rate, Q, of milk out of the alveoluscan then be found using Equation 3.6.Flow out of multiple alveoli.In order to determine the flow rate of milk from a collection of alveoli, through asystem of ducts, and out of the nipple, a similar process as for a single alveolus isused. To determine the flow rates in a branching network, Poiseuille’s law was appliedto each section of duct, with conservation of mass applied at junctions between twoor more pipes (flow into the junction must equal flow out of the junction). Eachalveolus is connected by a duct with radius rk and length lk to a main duct withradius rd,k and length ld,k between the kth junction and the next junction (k + 1 orthe nipple). Figure 2.2 shows the branching structure used and the location of eachparameter in the scheme.12Figure 2.2: Diagram of pressures and duct parameters in branching alveoli network.Black labels indicate prescribed parameters describing the branching network andgray labels indicate pressure values at a point in the network. The lower, horizontalduct is the main duct to the nipple (at the left end of the duct), while the verticalducts are branches between the alveoli and the main duct. Pbaby is the pressure atthe nipple from infant suckling.There are now n volume equations, similar to Equation 3.3, and n equations forconservation of mass and flow at each junction between ducts. This system of non-linear equations can then be used to solve for the alveoli radii, rA,k, and the pressureat the junctions, PJ,k, at each time step in order to determine the flow rate, Qk, outof each alveolus.Similarly to the single-alveolus case, the volume differential equation is quasi-steady-stated and each term can be rewritten as a function of rA,k. With multiplealveoli, the flow rate Qk will also depend on the junction pressure PJ,k. Following thesame derivation from Equation 3.3 to 3.7, we get fromdVA,kdt= −β (∆Pk −∆Πk)−Qk = 0 (3.8)to two cases: if rA,k ≥ r0, then Equation 3.8 becomes0 = −β KhrA,k(r2A,kr20− 1)+β(3Nk4pir3A,kRT − ΠC)− pir4k8νlk(KhrA,k(r2A,kr20− 1)+ Pω − PJ,k),(3.9a)13and if rA,k < r0, then Equation 3.8 becomes0 = β(3Nk4pir3A,kRT − ΠC)− pir4k8νlk(Pω − PJ,k). (3.9b)There are then n equations (3.9), one for each k = 1 to n. To solve for both rA,k andPJ,k, we need n more equations. In order to conserve mass in the mammary glandduct network, the net flux of milk at each junction of ducts must be zero. For eachalveolus, k = 2, . . . , n− 1, the net flow through the associated junction is0 =pir4k8νlk(Pk − PJ,k) +pir4d,k−18νld,k−1(PJ,k−1 − PJ,k)−pir4d,k8νld,k(PJ,k − PJ,k+1). (3.10a)For the first alveolus, there are only two ducts making up the junction (Figure 2.2),so for this junction, we have0 =pir418νl1(P1 − PJ,1)−pir4d,18νld,1(PJ,1 − PJ,2). (3.10b)Finally for the last alveolus, the duct with milk leaving the junction goes directly tothe infant, so the following pressure is PB, the pressure from infant suckling. For thisjunction, we then have0 =pir4n8νln(Pn − PJ,n) +pir4d,n−18νld,n−1(PJ,n−1 − PJ,n)−pir4d,n8νld,n(PJ,n − PB). (3.10c)The equations 3.9 and 3.10 can be solved simultaneously using a non-linear solver(fsolve in Matlab 2019a, Mathworks) at each time step to give the stretched radius,rA,k, and junction pressure, PJ,k, for each alveolus, k.2.3. Parameter ChoicesParameter values were chosen so that simulations produced physiologically rea-sonable results. Key characteristics that made a simulation “reasonable” include: allvariables remain non-negative (with the exception of pressure values which could benegative) and bounded, variables respond to initiation of feeding at an appropriatespeed (quickly for P , D and ω when not in a quasi steady state, and slowly for N ,G, C and r0), milk values decrease during feeding, C and/or r0 decrease in responseto increasing G, and fsolve succeeds to solve Equation 3.7 or Equations 3.9a-3.10cat each time-step. All parameter values used are given in the appendix in TablesA.1 and A.2. Unless otherwise stated, all simulations shown in the results use these14parameter values.153. Results3.1. Model 1: Constant milk removalNo milk age signalIn this model (Equations 1.1a-2.1a), aged milk is assumed to have no effect onmilk production and is thus not tracked separately from new milk. There is then onlyone slow variable: N , the amount of milk in the alveolus. Figure 3.1 shows a samplemilk profile for a single day with feeding once for one hour in every four hours.Figure 3.1: Simulated milk curve for 24 hours of feeding one hour in every four,as predicted by Model 1, constant milk removal with no milk aging. Gray regionsindicate time periods of feeding.As expected, the amount of milk oscillates: decreasing during feeding as theinfant removes milk and increasing between feedings as milk is produced. While thealveolus does not empty during feeding, it is reduced to less than half of the pre-feeding amount by the end of each feeding. Milk production has balanced so that theamount of milk produced matches the amount that is removed and the system hasno overall loss or gain of milk in the alveolus during this day.16Human breastmilk is about 7% lactose by weight, so if we use the molecularweight of lactose (342.3 g/mol) for the “dry” milk (N), this simulation correspondsto about 3.4 × 10−8 g of milk (including water weight) produced in a day (Jenness1979). Because only a single alveolus is being modeled and the same parameters willbe used for both constant removal and pressure-based removal in order to comparethe results more directly, this amount of milk produced is much smaller than that ofa mammary gland (Kent (2006) found that mothers produced 440 to 1220 g normallyin 24 hours). Differences between modeling alveoli of varying sizes will be discussedin a later section.To investigate the behavior of this model, the dynamics of the milk, N wereconsidered for two cases: between feedings (no milk removal) and during feeding(sustained milk removal). Figure 3.2 shows the phase portrait for N with this model.During a day with multiple feedings, oscillations in the amount of milk are causedby the system jumping between these two cases.17Figure 3.2: Phase lines for Model 1, constant milk removal with no milk aging, and thetwo system states: (a) between feedings (no milk removal) and (b) during feeding(sustained removal). Between feedings (a), there is a single stable steady state atN = 5.46×10−6 mol. With sustained feeding (b), there is a single stable steady stateat N = 2.97× 10−9 mol.Each case has a single stable steady state, with feeding moving the steady statefrom a relatively high value of 5.46×10−6 mol of milk down three orders of magnitudeto a relatively low 2.97 × 10−9 mol. The non-feeding steady state (a) dictates thelong-term behavior of the amount of milk after weaning. During the first stage ofinvolution, before structural changes to the mammary gland cause additional milkreabsorption, the system will move towards this value.Before weaning, the lactation system will oscillate in between these two steady18states, constantly jumping between the two phase lines. Notice that in the sampleday of feeding, milk remains between 0.5× 10−7 and 2× 10−7 mol—near neither ofthe two steady states (Figure 3.2).Although the during feeding equilibrium gives the long-term result of milk in themammary gland with no feeding, this model can only describe up to two days afterfeeding ends—through the first stage of involution when production slows. After twodays, structural changes occur and leftover milk is reabsorbed, guided by hormonaleffects not described in this model (Jindal et al. 2014; Watson 2006). . Figure 3.3shows milk production during the two days after feeding ends.Figure 3.3: Simulated milk curve for the end of feeding, as predicted by Model 1,constant milk removal without milk aging. Periods of feeding are shown in gray andend at hour 24.After feeding ends, milk rises quickly after the final feeding and is ten timeshigher by the end of the two days of no feeding. However, milk does not reach theequilibrium point of 5.46 × 10−6 mol during this time frame. In fact, it does notapproach this equilibrium until about four weeks after feeding ends, far beyond thescope of this model.Milk age affecting cell quiescenceIn this model, new milk ages and is converted to aged milk, G, and causesalveolar cell quiescence (Equations 1.1a-2.3). There are now three slow variables: N ,19new milk, G, aged milk, and C, proportion of alveolar cells active. Figure 3.4 showsa sample day with feeding occurring regularly once in every four hours.Figure 3.4: Simulated milk curve for 24 hours of feeding one hour in every four, aspredicted by Model 1, constant milk removal with aged milk affecting cell quiescence.Gray regions indicate time periods of feeding.New and aged milk both oscillate, increasing between feedings and decreasingduring feeding. The fraction of cells active, however, increases while aged milk is lowand decreases when aged milk is sufficiently high, causing oscillations with the sameperiod as for milk. Over all, aged milk is relatively low and the fraction of cells activeremains close to one.Including the effect of milk aging on cell quiescence did not change the dailypattern of milk production during feeding, but it may change milk production duringinvolution. Figure 3.5 shows the lactation system in the two days after the finalfeeding compared to the last day of feeding.20Figure 3.5: Simulated milk curve for the end of feeding, as predicted by Model 1,constant milk removal with aged milk affecting cell quiescence. Periods of feeding areshown in gray and end at hour 24.In this case, production of new milk slows soon after the final feeding while agedmilk production speeds up. By the end of the two days, new milk has leveled off whileaged milk continues to grow quickly. The fraction of cells active drops as aged milkrises and continues to decrease steadily to the end of the two day period. AlthoughC decreases slowly and does go below 85% of cells active, this small decrease issufficient to slow new milk production to matching the rate of aging within the twoday involution period.Milk age affecting alveolus sizeIn this model, aged milk causes alveolar cell death (Equations 1.1a-1.5b, 2.1b,2.2, 2.4). The slow variables are then N , G and r0, the reference (non-stretched)radius of the alveolus. Figure 3.6 shows a sample day of regular feeding with thismodel.21Figure 3.6: Simulated milk curve for 24 hours of feeding one hour in every four, withModel 1, constant milk removal with aged milk affecting alveolus size. Gray regionsindicate time periods of feeding.New and aged milk oscillate similarly to the previous model. The referenceradius also oscillates in a similar pattern to C in the previous model—increasingwhile aged milk is low and decreasing once aged milk is sufficiently high.Involution also shows similarities between the two models. Figure 3.7 shows thelactation system during involution.22Figure 3.7: Simulated milk curve for the end of feeding, as predicted by Model 1,constant milk removal with aged milk affecting alveolus size. Periods of feeding areshown in gray and end at hour 24.The alveolus radius decreases later and more gradually after feeding ends thanthe fraction of cells active decreased. However, the same effect on new and aged milkproduction is seen: new milk production slows until new milk levels off while agedmilk rises rapidly until the end of the period. Additionally, although r0 decreasesmore gradually than C, in this case new milk begins to decrease slightly by the endof involution stage one.Milk age affecting cell quiescence and alveolus sizeIn this final version of Model 1, aged milk both deactivates alveolar cells andcauses alveolar cell death. All four slow variables are in effect: N , G, C and r0. Figure3.8 shows a sample day of milk production with regular feeding occurring for onehour in every four hours.23Figure 3.8: Simulated milk curve for 24 hours of feeding one hour in every four, withModel 1, constant milk removal with aged milk affecting cell quiescence and alveolussize. Gray regions indicate time periods of feeding.When paired together, the effects of milk age on active cells and alveolus radiusmaintain the same patterns as seen when they were modeled separately. The fractionof cells active still oscillates between 0.995 and 0.999, increasing while aged milk is lowand decreasing when it is high (Figure 3.8). Similarly, the reference radius oscillatesnear 9.97× 10−3 mm both when it is included alone and with C (Figure 3.8). Theseoscillations follow the aged milk pattern, which is also maintained. New and agedmilk increase between feedings and decrease during feeding, with aged milk oscillatingbetween 1× 10−8 mol and 3× 10−8 mol, and new milk oscillating between 5× 10−8mol and 1.5× 10−7 mol (Figure 3.8). Involution also shows similarities between thetwo models. Figure 3.9 shows the lactation system during involution.24Figure 3.9: Simulated milk curve for the end of feeding, as predicted by Model 1,constant milk removal with aged milk affecting cell quiescence and alveolus size.Periods of feeding are shown in gray and end at hour 24.Again, the same general patterns are found when the two milk aging effectsare paired together. All four versions of Model 1 show similar daily patterns and allshow marked increases in milk levels during the first stage of involution, althoughnot including the effects milk aging causes milk to rise fastest after the final feeding.Including the effect of milk aging in any capacity produced similar dynamics in milkproduction as well as the variables affected by milk age.Cluster feeding with constant milk removalIncreased frequency of feeding has been shown to increase production in womenwith low milk supply, so various cluster feeding schedules were tested on each version25of Model 1 and the amount of milk produced compared (Fewtrell et al. 2016; Dalyet al. 1996). The previously used feeding schedule of feeding one hour in every fourwas then modified to test different numbers of feedings in a single cluster. The amountof time spent feeding per period (one hour) and the length of the period (four hours)were maintained, but the hour of feeding was split into Nc feedings 20 minutes apart.Figure 3.10 shows the amount of milk produced for the four constant milk removalmodels discussed in this section with six different feeding schedules (Nc = 1 is theregular feeding schedule used previously, Nc > 1 are cluster feeding schedules asdescribed). The amount of milk produced was measured by summing the productioncomponent of the N equation (Equation 2.1b) across all time steps.Figure 3.10: Total milk produced in 24 hours for each of the four versions of Model 1,and six different feeding schedules. A single feeding in a cluster (Nc = 1) correspondsto a regular feeding schedule with one hour of feeding in every four hours. Highernumbers of feedings in a cluster correspond to one hour of feeding split into Ncfeedings, 20 minutes apart, every four hours.26For each model, the amount of milk produced increases with the number offeedings in a cluster, although only slightly. Both versions of the model in which milkaging affects alveolus size (r0) produce consistently higher milk across all numbers offeeding. Including the effect of milk aging on the proportion of active cells (C) lowersthe amount of milk produced. In cases where age affects r0, the alveolus is allowedto grow in response to low levels of aged milk, generating a higher milk productioncapacity. In the cases where age affects C, the alveolar cells are allowed to deactivatein response to high levels of aged milk rather than having C fixed at one. Eventhough including the effect of aged milk on C decreased production, it also increasedthe system’s sensitivity to changes in the feeding schedule, as seen by the differencesin slope in the bottom two curves.Although continuing to increase the number of feedings per cluster at somepoint generates an entirely impractical feeding schedule, methods of increasing milkproduction are of interest. Thus the time in between feedings in a cluster, previouslyfixed at 20 minutes, was varied in order to allow for more feedings in a cluster whilekeeping the period of the feeding schedule fixed at four hours. Figure 3.11 shows theamount of milk produced in 24 hours for all four models as the time between clustersand number of feedings in a cluster were varied.27Figure 3.11: Total milk produced in 24 hours for each of the four versions of Model 1,with varied feeding schedules parameterized by the time between feedings in a singlecluster (tb) and the number of feedings in a cluster (Nc). The feeding schedule periodwas held fixed at 4 hours and the total amount of time spent feeding per period wasfixed at one hour.Each model shows a similar pattern: the time between feedings in cluster (tb)does not noticeably affect the amount of milk produced while the number of feedingsin a cluster (Nc) does change the amount of milk produced slightly. Although thevariations in milk produced for each Nc seems random, the pattern persists acrossall models. Although this persistent, “random” pattern may be the result of coarsetime-steps generating errors in approximating the feeding schedule, use of a smallertime step produced a similar pattern. Figure 3.12 shows the milk production for thefinal Model 1 version (aged milk affecting C and r0) using a time-step of 0.01 min.28Figure 3.12: Total milk produced in 24 hours as predicted by Model 1 with agedmilk affecting cell quiescence and alveolus size, and with varied feeding schedulesparameterized by the time between feedings in a single cluster (tb) and the numberof feedings in a cluster (Nc). A smaller time-step of 0.01 min was used to reduce errorin the feeding schedule approximation. The feeding schedule period was held fixedat 4 hours and the total amount of time spent feeding per period was fixed at onehour.Refining the time-step smooths the pattern somewhat and distinguishes betweenthe amount of milk produced between different times between feedings in a cluster.However, there are still drops in the amount of milk produced as Nc increases. It maybe that the difference in milk production between one feeding schedule and the next issmall enough that even very small errors in the feeding schedule approximation couldlead to significant differences in the amount of milk produced. This could explain thepersistent peaks at Nc = 10, 12, 15, 20 and 30, which all evenly divide 60 minutes.Due to substantial computation times, time steps smaller than 0.01 min were not29explored for this section. In order to determine what might cause the decreases inproduction, two similar feeding schedules were compared: Nc = 20 and Nc = 21,with tb = 5 min used for both. Figure 3.13 shows sample days of milk productionwith these two feeding schedules.Figure 3.13: Simulated milk curve for 24 hours with Model 1, constant milk removalwith aged milk affecting cell quiescence and alveolus size. Clusters of feeding occurevery four hours with 5 minutes between feedings in a cluster, one hour of feedingtotal for each cluster and either (a) 20 or (b) 21 feedings in a cluster. Gray regionsindicate time periods of feeding. With (a) 20 feedings in a cluster 8.2367 × 10−6mol of milk is produced in these 24 hours, while with (b) 21 feedings in a cluster8.1640× 10−6 mol is produced.With Nc = 21, both new and aged milk levels are consistently higher than withNc = 20, which in turn causes both C and r0 to be lower. Thus milk productionis lower with 21 feedings in a cluster. The amount of time spent feeding is likelyunderestimated then—because new and aged milk levels are higher but the amount ofmilk produced is lower, less milk must be removed during feeding than with Nc = 20.With less milk removed, less milk is produced and the lactation system tends towardan equilibrium.Overall, there is a general increase in milk production as the number of feedingsin a cluster increases and as the time between feedings in a cluster increases. Taking30only values of Nc with 60 mod Nc = 0, milk production strictly increases with thenumber of feedings per cluster. Increasing either Nc or tb moves the feeding sched-ule closer to a regular feeding schedule (all feedings equally spaced) with a higherfrequency of feeding than the original regular feeding schedule used.3.2. Model 2: Pressure-based milk removalIn order to compare modeling a single alveolus to multiple alveoli, the flow rateof milk out of the alveoli must be based on the physics and geometry of the mammarygland rather than simply having a constant, prescribed flow rate. Section 2.2 describeshow this flow rate is calculated for each case: a single alveolus and multiple alveoli.Because of the requirement of solving a non-linear equation or system of equationsat each time-step, times for completely simulations are fairly slow, taking about 40minutes to simulate 36 hours of lactation. For both cases, the final version of Model1, constant milk removal with aged milk affecting both cell quiescence and alveolussize, was used to model milk production.Flow out of a single alveolusBefore testing multiple alveoli, the behavior of the model with the new pressure-based milk removal was investigated. Using the same parameters from Model 1, aswell as new, physical parameters for the fluid flow, one day of feeding regularly for onehour in every four was simulated. Figure 3.14 shows the simulated milk productionfor a day early in lactation onset while the alveolus is still growing.31Figure 3.14: Simulated milk curve for 24 hours of feeding one hour in every four, withModel 2, pressure-based milk removal with a single alveolus. Gray regions indicatetime periods of feeding.Although new and aged milk oscillate similarly with constant or pressure-basedremoval, pressure-based removal necessarily led to emptying of the alveolus at eachfeeding. When parameters that may lead to slower milk removal were used, either thenumeric solver (Matlab’s fsolve) failed to solve Equation 3.7 (e.g., larger duct radiusrd or larger maximum reference radius r0,max) or large changes in the parameter hadlittle effect and did not prevent the alveolus from emptying (e.g., changing elasticmodulus K of the alveolus wall or the pressure PB of infant suckling). In this case,pressure build up inside the alveolus as well as oxytocin-induced contraction of thealveolus wall lead to immediate release of milk upon suckling, demonstrating thelet-down reflex (Newton and Newton 1948).The fraction of cells active oscillate near one, as in Model 1, with the sameperiod. However, with aged milk low and removed immediately upon feeding, C is32consistently higher and experiences much sharper troughs in its oscillation. There isno smooth increase in C after feeding begins and milk is removed; instead C imme-diately increases upon the beginning of feeding and C spends more time increasingthan decreasing. Although the alveolus reference radius is still growing, milk levelsare so low due to the regular emptying of the alveolus that milk production is notsignificantly affected by the lowered reference radius. Because involution is charac-terized by the lack of feeding and because implementation of feeding is the differencebetween Models 1 and 2, involution will follow the same pattern in this model as inthe final version of Model 1 and is not shown.Different cluster feeding schedules were applied to this model and amounts ofmilk produced compared. To reduce computation time, the time-step was not refinedfurther. Figure 3.15 shows the amount of milk produced in 24 hours for each clusterfeeding schedule.33Figure 3.15: Total milk produced in 24 hours for each of the four versions of Model 2in a single alveolus, with varied feeding schedules parameterized by the time betweenfeedings in a single cluster (tb) and the number of feedings in a cluster (Nc). Thefeeding schedule period was held fixed at 4 hours and the total amount of time spentfeeding per period was fixed at one hour.As in Model 1, milk production generally increases with more feedings in acluster, with some sporadic decreases possibly caused by errors in approximatingthe feeding schedule. In this case, the time between feedings in a cluster has aneffect even for the coarse time-step, with milk production increasing slightly as tbincreases. Increasing either Nc or tb moves the feeding schedule closer to a regularfeeding schedule (all feedings equally spaced) with a higher frequency of feeding thanthe original regular feeding schedule used.Including pressure-based milk removal necessitates alveolus emptying and gen-34erates a system very sensitive to parameters describing the geometry of the simulatedmammary gland. Additionally, increases in milk production with increasing time be-tween feedings in a cluster were present, and therefore greater milk production inmore regular and frequent feeding schedules.Flow out of multiple alveoliAlthough the number of alveoli in a mammary gland is far larger than three(in dogs, there can be as many as 3,000 alveoli), to save on computation time, threealveoli were modeled (Orfanou et al. 2010). The general geometry of the three alveoliwas as shown in Figure 2.2. Each alveolus was 1 mm away from the main duct (lk = 1mm for k = 1, 2, 3) and the main duct was 6 mm long, with alveoli evenly spacedalong the duct so that ld,k = 2 mm for k = 1, 2, 3. All duct radii were equal withrk = rd,k = 0.0005 mm for k = 1, 2, 3 (larger duct radii caused fsolve to fail to solveEquations 3.9a-3.10c for the junction pressures PJ,k and stretched alveolus radii rA,k).35Figure 3.16: Simulated milk curve for 24 hours of feeding one hour in every four, withModel 2, pressure-based milk removal with three alveoli. Alveolus 1 is furthest fromthe nipple while alveolus 3 is closest. Gray regions indicate time periods of feeding.Each of the three alveoli behave the exact same way, and behave the same asin the single alveolus case (Figure 3.14). Additional milk volume entering the mainduct has not slowed milk removal or prevented the alveoli from emptying immediatelyupon feeding onset. Comparing the application of different feeding schedules showssimilar patterns as well. Figure 3.17 shows the amount of milk produced in totalfor all three alveoli for different feeding schedules and a reduced number of samplevalues for Nc to save computation time.36Figure 3.17: Total milk produced in 24 hours for each of the four versions of Model2 in three alveoli, with varied feeding schedules parameterized by the time betweenfeedings in a single cluster (tb) and the number of feedings in a cluster (Nc). Thefeeding schedule period was held fixed at 4 hours and the total amount of time spentfeeding per period was fixed at one hour.With three alveoli instead of one, the total capacity for amount of milk pro-duced is tripled and this is reflected in the amount of milk produced for each feedingschedule, which is about triple the amount produced in the single alveolus case (Fig-ure 3.17). As with the single alveolus cases, the amount of milk produced generallyincreases as feedings become more frequent and regular (Figures 3.11, 3.15).Although the dynamics of milk production are not noticeably different whenmodeling more alveoli, the amount of milk produced is greater with multiple alveoli.In Figure 3.16 the total amount of milk produced in the displayed 24 hours is 2.455×10−6 mol. To compare, a single alveolus was simulated with a maximum volume equalto the total maximum volume of three alveoli (r0,max = 0.01 · 31/3 mm for the single37alveolus and r0,max = 0.01 mm for the three alveoli). The single alveolus produced1.706× 10−6 mol in 24 hours, about 2/3 of the amount produced with three alveoli.Because r20 directly affects the milk production rate (Equation 2.1b), but also affectsthe amount of stretch experienced by an alveolus, as well as the pressure Pk insidethe alveolus, the amount of milk produced is dependent upon alveolus surface area aswell as alveolus volume. Thus modeling more alveoli allows for greater surface areaand greater milk production.Because the alveoli were emptied immediately upon feeding in both the singlealveolus and multiple alveolus cases, each of the three alveoli are exact copies of eachother and demonstrate no change in the dynamics of milk production. In order to seea noticeable difference between alveoli, the distance between alveoli along the mainduct needs to be large: ld,k = 30 mm. Figure 3.18 shows a day of regular feeding withthree alveoli each 30 mm apart.38Figure 3.18: Simulated milk curve for 24 hours of feeding one hour in every four,with Model 2, pressure-based milk removal with three alveoli. Alveolus 1 is furthestfrom the nipple while alveolus 3 is closest, with ld = 30 mm between each alveolusalong the main duct. Gray regions indicate time periods of feeding.Increasing the distance between alveoli to 30 mm causes the alveolus closestto the nipple (alveolus 3, shown in yellow) to empty faster. Alveolus 1 (blue), thefurthest from the nipple, continues to gain milk at the onset of feeding until boththe closer alveoli are emptied. The fractions of cells active and alveoli reference radiiare now noticeably different between alveoli, but only very slightly. As milk is beingproduced between feedings, the alveoli are producing milk seemingly identically, andit is only once feeding begins that there alveoli become distinct, suggesting that thesmall changes to C and r0 between alveoli are not large enough to significantly affectmilk production. It is possible that even longer, less physiological, duct lengths could39allow for greater differentiation between alveoli, but this would necessitate folding ofthe duct in order for it to fit inside the mammary gland and imaging of the breastdoes not show this kind of structure (Ramsay et al. 2005).404. Discussion4.1. Model 1: Constant milk removalWith Model 1, two different local effects for aged milk were explored: agedmilk deactivating alveolar cells and aged milk causing alveolar cell death. These twoeffects, together and separately, did not noticeably change the daily milk patternsand caused only minor changes to the pattern of involution. Although including theeffects of aged milk on production did cause new milk to level off fairly quickly afterfeeding ended in the first stage of involution, aged milk continued to rise rapidly.Thus the total amount of milk, new and aged, still rose rapidly throughout thefirst stage of involution. However, the marked loss of productive cells due to milkaging, either by quiescence or apoptosis of alveolar cells, is a key marker of the firststage of involution (Watson 2006; Geddes 2007; Schere-Levy et al. 2003). While thedifferent aging effects did not necessarily change how total milk levels behaved duringinvolution, both milk aging effects captured a key hallmark of involution stage one(Figures 3.1-3.9).The main difference between the aging effects was the total amount of milkproduced in 24 hours (Figures 3.10-3.12). Because when the fraction of alveolar cellsactive C was not changing with aged milk C was held constant at one, allowingC to change could only decrease the alveolus’ capacity for milk production. Thusincluding the effect of aged milk on the fraction of cells active always decreased thetotal amount of milk produced. The alveolus reference radius, however, was set belowits maximum when not changing with aged milk. Thus allowing the reference radiusto change meant the alveolus could grow as well as shrink, and the capacity for milkproduction could increase. Thus holding both C and r0 at their maximums when theywere not changing would mean that both could cause decreases in milk productionwhen allowed to change with aged milk.These particular values were chosen to best match the cases in which eachvariable would not be affected by aged milk. If the fraction of cells active is not41affected by aged milk and quiescent cell activation is not a method of modulatingmilk production, then the presence of quiescent/stem cells likely do not influence theproduction capacity of a single alveolus, although they may cells may be activatedor deactivated in order to change the total number of alveoli (Hassiotou and Geddes2013; Orfanou et al. 2010). For simplicity, it was then assumed that in these cases thefraction of active cells was one. If the reference radius of the alveolus is not affectedby aged milk, then there must be some other method of changing alveolus size be-cause alveoli can change size throughout lactation (Orfanou et al. 2010). Because weare interested in the low milk production case and methods of increasing production,it was assumed that alveoli were not at their maximum size. Using different assump-tions or considering different cases would allow milk production to either increase ordecrease with the addition of aged milk effects. Further research into the dynamics ofmammary cells during lactation is needed to determine which assumptions are mostrelevant when considering mothers with low production.Milk aging could affect milk production through either of the considered mech-anisms. Both were capable of affecting milk production and both exhibited similardynamics. Additionally both mechanisms are equally capable of explaining how in-sufficient breast drainage and prolonged periods between feedings could decreasemilk production. If milk aging causes mammary cells to deactivate, then accumula-tion of aged milk—from either insufficient milk drainage during feeding or prolongedabsence of feeding—would cause alveolar cells to become quiescent and stop pro-ducing milk until enough aged milk could be removed to allow quiescent cells tobe re-activated. If milk aging causes mammary cell apoptosis, then accumulation ofaged milk would cause milk-producing cells to die, thereby decreasing the capacityfor production until enough aged milk could be removed to decrease the cell deathrate and allow the cell population to increase. Because both aged milk effects werefound to be biologically feasible and to reproduce known lactation dynamics equallywell, both effects were used in Model 2.424.2. Model 2: Pressure-based milk removalIncluding fluid dynamics and pressure-based milk removal into feeding withModel 2 allowed for the exploration of the effects of the geometry of the mammarygland on milk production. This model, however, was very sensitive to certain pa-rameters, in particular the radius of the ducts rd,k, and seemed to have little to nodependence on other parameters, most notably the pressure of infant suckling PB.Under this restrictive parameter regime, the alveoli always empty during feed-ing (Figures 3.14, 3.16, 3.18). This forced emptying generated sharper troughs inthe oscillations of the fraction of cells active, and prevented any differences in milkproduction from emerging. When multiple alveoli were simulated, they were thenmore like copies of the single alveolus case rather than a separate model to consider.Extending the duct length beyond what is physiologically realistic may allow for newpatterns to emerge. Additionally, modeling more alveoli, allowing for a more complexgeometry, may reveal differences in milk production between alveoli and would be amore realistic model to use for investigating low milk production.The mammary gland could be approximated as a single large alveolus ratherthan many small alveoli, and this approximation was used when initially creatingModel 1. However, collecting the volume of all the alveoli in the breast together intoa single larger alveolus reduces milk production capacity. As discussed in section 3.2,milk production depends upon alveolus surface area as well as volume, so increasingthe total amount of surface area in the mammary gland allows for higher milk pro-duction. This may be why the mammary gland evolved a structure more similar tothe lungs, with many small alveoli, than the stomach, with one large sac. One largesac is energetically less expensive to develop—fewer cells need to be made and keptup, and it is easier to deliver nutrients to each cell with blood vessels. The fact thatmammals have evolved branched mammary gland structures with many alveoli sug-gests that high milk production in mothers is essential to the success of mammalianoffspring.434.3. Increasing milk productionBecause low milk production has become an issue for many women, methodsof increasing production must be found and changing the feeding schedule has beenproposed as a solution (Kent et al. 2012; Daly et al. 1996). Using the models pre-sented, we were able to determine that increased frequency of feeding does increasemilk production. Additionally, more time feeding also increases the amount of milkproduced. Splitting a single feeding into a cluster of multiple feedings did increaseproduction, but production increased further if those feedings were spaced out and nolonger clustered. Therefore cluster feeding can cause an increase in milk production,but only through the increase in the number of feedings in a day. The optimal feed-ing schedule seems to be feeding as frequently, as regularly and as much as possible(Figures 3.11, 3.12, 3.15 and 3.17).Unfortunately, the most frequent feeding schedules examined in Figures 3.11,3.12, 3.15 and 3.17 are incredibly unrealistic for any human being. Feeding for twominutes every five minutes is impractical and does not allow for much time to doinganything other than feeding—if the infant can even be convinced to feed in thosetwo minutes. This “optimal” feeding schedule is then not really optimal and couldnot reasonably be suggested as a method for increasing milk production. In orderto determine and truly optimal feeding schedule, many other components of infantgrowth and health would need to be included, such as sleep. Cluster feeding may bea way to increase the number of feedings in a day while still allowing time for bothmother and infant to do other things.4.4. Future workIn order to fully explore the presented models, more alveoli could to be simu-lated. This model could also be combined with a model similar to the fluid dynamicsmodel of the breast by Negin Mortazavi et al. (2016). With a more complete simula-tion of the actual mammary gland, new milk production patterns may emerge, andthe feeding schedules tested may have a different effect on production.This model could also be paired with existing metabolism models, like the modelby Hall (2006), to track infant growth throughout lactation. Such a model could44explore the effects of hunger and compare infant growth between the infant demandfeeding schedule and more regular feeding schedules. If combined with circadianrhythm and sleep models, the development of a normal (adult-like) sleep cycle couldbe studied.Although the reasons why low milk production is common in developed coun-tries is unknown, factors like stress are known to affect milk production (Dewey2001). Other biological signalling molecules are also known to affect lactation, in-cluding TRH, serotonin, insulin and insulin-like growth factor (IGF) (Lemay et al.2013; van de Moosdijk et al. 2017; Horseman and Collier 2013). Model 1 could bemodified to include these factors in order to investigate possible causes of low milkproduction, and therefore possible pharmacological interventions.Because not all physical and biological components of lactation were included,this model cannot describe lactation as a whole. However, for investigating the ef-fects of different feeding schedules on milk production, the presented models balancesimplicity with predictive power and greater complexity is not necessary.45ReferencesAmis, D. 2015. A Childbirth Educators Commentary on Hormonal Physiology ofChildbearing: Evidence and Implications for Women, Babies, and Maternity Care. The Journal of Perinatal Education, 24(3):154–159. [Cited on page 1.]Bertram, R., Egli, M., Toporikova, N., and Freeman, M. E. 2006. A mathematicalmodel for the mating-induced prolactin rhythm of female rats. American Journalof Physiology-Endocrinology and Metabolism, 290(3):573–582. [Cited on page 5.]Buckley, S. J. 2015. Executive Summary of Hormonal Physiology of Childbearing:Evidence and Implications for Women, Babies, and Maternity Care. The Journalof Perinatal Education, 24(3):145–153. [Cited on page 1.]Cox, D. B., Owens, R. A., and Hartmann, P. E. 1996. Blood and milk prolactin andthe rate of milk synthesis in women. Experimental Physiology, 81(6):1007–1020.[Cited on pages 3 and 4.]Dag, B., Keskin, I., and Mikailsoy, F. 2005. Application of different models to thelactation curves of unimproved Awassi ewes in Turkey. South African Journal ofAnimal Sciences, 35(4):238–243. [Cited on page 2.]Daly, S., Owens, R., and Hartmann, P. 1993. The shortterm synthesis and infantregu-lated removal of milk in lactating women. Experimental Physiology, 78(2):209–220.[Cited on page 3.]Daly, S. E., Kent, J. C., Owens, R. A., and Hartmann, P. E. 1996. Frequencyand degree of milk removal and the short-term control of human milk synthesis.Experimental Physiology, 81(5):861–875. [Cited on pages 1, 3, 26, and 44.]Dewey, K. G. 2001. Maternal and fetal stress are associated with impaired lactoge-nesis in humans. The Journal of nutrition, 131(11):3012S–5S. [Cited on pages 1and 45.]46Dijkstra, J., France, J., Dhanoa, M., Maas, J., Hanigan, M., Rook, A., and Beever,D. 2010. A Model to Describe Growth Patterns of the Mammary Gland DuringPregnancy and Lactation. Journal of Dairy Science, 80(10):2340–2354. [Cited onpages 2 and 4.]Egli, M., Leeners, B., and Kruger, T. H. 2010. Prolactin secretion patterns: Basicmechanisms and clinical implications for reproduction. Reproduction, 140(5):643–654. [Cited on page 2.]Ferreira, A. G., Henrique, D. S., Vieira, R. A., Maeda, E. M., and Valotto, A. A.2015. Fitting mathematical models to lactation curves from holstein cows in thesouthwestern region of the state of Parana, Brazil. Anais da Academia Brasileirade Ciencias, 87(1):503–518. [Cited on page 2.]Fewtrell, M. S., Kennedy, K., Ahluwalia, J. S., Nicholl, R., Lucas, A., and Burton,P. 2016. Predictors of expressed breast milk volume in mothers expressing milkfor their preterm infant. Archives of Disease in Childhood: Fetal and NeonatalEdition, 101(6). [Cited on pages 1 and 26.]Geddes, D. T. 2007. Inside the Lactating Breast: The Latest Anatomy Research.Journal of Midwifery and Women’s Health, 52(6):556–563. [Cited on page 41.]Grattan, D. R. 2015. 60 YEARS OF NEUROENDOCRINOLOGY: Thehypothalamo-prolactin axis. Journal of Endocrinology, 226(2):T101–T122. [Citedon page 2.]Hall, K. D. 2006. Computational model of in vivo human energy metabolism duringsemistarvation and refeeding. American Journal of Physiology-Endocrinology andMetabolism, 291(1):E23–E37. [Cited on page 44.]Hassiotou, F. and Geddes, D. 2013. Anatomy of the human mammary gland: Currentstatus of knowledge. Clinical Anatomy, 26(1):29–48. [Cited on page 42.]47Horseman, N. D. and Collier, R. J. 2013. Serotonin: A Local Regulator in theMammary Gland Epithelium. Annual Review of Animal Biosciences, 2(1):353–374. [Cited on pages 8, 9, and 45.]Hurst, N. M. 2007. Recognizing and Treating Delayed or Failed Lactogenesis II.Journal of Midwifery and Women’s Health, 52(6):588–594. [Cited on page 1.]Jenness, R. 1979. The composition of human milk. Seminars in perinatology,3(3):225–239. [Cited on page 17.]Jiang, H. and Sun, S. X. 2013. Cellular pressure and volume regulation and impli-cations for cell mechanics. Biophysical Journal, 105(3):609–619. [Cited on page11.]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 involu-tion reveals regression of secretory lobules mediated by tissue-remodeling. BreastCancer Research, 16(2):1–14. [Cited on page 19.]Kent, J. C. 2006. Volume and Frequency of Breastfeedings and Fat Content of BreastMilk Throughout the Day. Pediatrics, 117(3):e387–e395. [Cited on page 17.]Kent, J. C., Prime, D. K., and Garbin, C. P. 2012. Principles for Maintaining or In-creasing Breast Milk Production. Journal of Obstetric, Gynecologic, and NeonatalNursing, 41(1):114–121. [Cited on pages 1, 3, and 44.]Lemay, D. G., Ballard, O. A., Hughes, M. A., Morrow, A. L., Horseman, N. D., andNommsen-Rivers, L. A. 2013. RNA Sequencing of the Human Milk Fat Layer Tran-scriptome Reveals Distinct Gene Expression Profiles at Three Stages of Lactation.PLoS ONE, 8(7). [Cited on page 45.]Marasco, L. A. 2015. Unsolved Mysteries of the Human Mammary Gland: Definingand Redefining the Critical Questions from the Lactation Consultants Perspective.Journal of Mammary Gland Biology and Neoplasia, 19(3-4):271–288. [Cited onpage 1.]48Mathworks, C. 2019. Function Reference R 2019 a. [Cited on pages 12 and 14.]Negin Mortazavi, S., Geddes, D., and Hassanipour, F. 2016. Lactation in the Hu-man Breast From a Fluid Dynamics Point of View. Journal of BiomechanicalEngineering, 139(1):011009. [Cited on pages 2 and 44.]Newton, M. and Newton, N. R. 1948. The let-down reflex in human lactation. TheJournal of Pediatrics. [Cited on pages 11 and 32.]Orfanou, D. C., Pourlis, A., Ververidis, H. N., Mavrogianni, V. S., Taitzoglou, I. A.,and Boscos, C. M. 2010. Histological Features in the Mammary Glands of FemaleDogs throughout Lactation. pages 473–478. [Cited on pages 35 and 42.]Ramsay, D. T., Kent, J. C., Hartmann, R. A., and Hartmann, P. E. 2005. Anatomyof the lactating human breast redefined with ultrasound imaging. Journal ofAnatomy, 206(6):525–534. [Cited on page 40.]Rattanakul, C. and Lenbury, Y. 2009. A mathematical model of prolactin secre-tion: Effects of dopamine and thyrotropin-releasing hormone. Mathematical andComputer Modelling, 49(9-10):1883–1892. [Cited on page 5.]Schere-Levy, C., Buggiano, V., Quaglino, A., Gattelli, A., Cirio, M. C., Piazzon, I.,Vanzulli, S., and Kordon, E. C. 2003. Leukemia inhibitory factor induces apop-tosis of the mammary epithelial cells and participates in mouse mammary glandinvolution. Experimental Cell Research. [Cited on pages 3, 8, and 41.]Schmitt-Ney, M., Happ, B., Hofer, P., Hynes, N. E., and Groner, B. 2014. Mammarygland-specific nuclear factor activity is positively regulated by lactogenic hormonesand negatively by milk stasis. Molecular Endocrinology, 6(12):1988–1997. [Citedon pages 3 and 8.]Sunaric´, S., Jovanovic´, T., Spasic´, A., Denic´, M., and Kocic´, G. 2016. ComparativeAnalysis of the Physicochemical Parameters of Breast Milk, Starter Infant For-mulas and Commercial Cow Milks in Serbia. Acta Facultatis Medicae Naissensis,33(2):101–108. [Cited on pages 10 and 51.]49Tay, C. C., Glasier, A. F., and McNeilly, A. S. 1996. Twenty-four hour patternsof prolactin secretion during lactation and the relationship to suckling and theresumption of fertility in breast-feeding women. Human Reproduction, 11(5):950–955. [Cited on pages 2, 6, and 7.]van de Moosdijk, A., Fu, N., Rios, A., Visvader, J., and Van Amerongen, R. 2017.Lineage Tracing of Mammary Stem and Progenitor Cells., volume 1501. [Citedon page 45.]Walker, S. P. 1997. Nutritional issues for women in developing countries. Proceedingsof the Nutrition Society, 56(1B):345–356. [Cited on page 1.]Watson, C. J. 2006. Key stages in mammary gland development Involution: Apop-tosis and tissue remodelling that convert the mammary gland from milk factoryto a quiescent organ. Breast Cancer Research, 8(2):1–5. [Cited on pages 3, 19,and 41.]Wilde, C. J. and Peaker, M. 1990. Autocrine control in milk secretion. The Journalof Agricultural Science, 114(3):235–238. [Cited on pages 3, 6, and 8.]Wood, P. 1968. Factors affecting persistency of lactation in cattle. Nature, 218:894.[Cited on page 2.]50A. Parameter ValuesParameter Value used DescriptionR 8.31441 J/K/mol Ideal gas constantT 310 K Body temperatureν 2.4× 10−3 MPa·s Viscosity of milk (Sunaric´ et al. 2016)Table A.1: Milk production model physical constants. These values were not changedthroughout the model development process.Parameter Value used Parameter Value usedkP1 1 min−1 kP2 0.5kP3 2 min−1 kD1 0.25 min−1kD2 0.75 kD3 0.5kD4 0.66 min−1 kD5 2.5 min−1kU 1 kS1 1 MPa−1kS2 16.5 MPa kσ 1 MPa/molkω1 1 min−1 kω2 1 min−1kN 10−5 mol/mm2/min kl 0.0001 min−1kFl 0.2 min−1 kG 0.001 min−1kC1 0.01 min−1 kC2 103 mol−1min−1kr1 2 mm−2min−1 kr2 1 mm−2min−1r0,max ∗ 0.01 mm kr3 700 mol−1min−1ld,ld,k,lk ∗ 1 mm rd,rd,k,rk 0.0005 mmK 100 MPa h 0.0005 mmkPω 2× 10−9 MPa PB (Pbaby) 0.001 MPaβ 1 ΠC 0.000625 mol/mm3Nc ∗ 1 tb ∗ 20 min∆t ∗ 1 minTable A.2: Milk production model parameters. See Model chapter for descriptionof parameters. These parameters were modified throughout the model developmentprocess and final parameter values were chosen based on whether simulations runwith these parameters produced physiologically reasonable results. Parameters witha star (∗) indicate those that were changed for tests shown in this thesis. Furtherdiscussion of parameter choices is included in the Model, Results and Discussionchapters.51


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items