A Computational Model of the Rainbow Trout Hypothalamus-Pituitary-Ovary-Liver Axis

Reproduction in fishes and other vertebrates represents the timely coordination of many endocrine factors that culminate in the production of mature, viable gametes. In recent years there has been rapid growth in understanding fish reproductive biology, which has been motivated in part by recognition of the potential effects that climate change, habitat destruction and contaminant exposure can have on natural and cultured fish populations. New approaches to understanding the impacts of these stressors are being developed that require a systems biology approach with more biologically accurate and detailed mathematical models. We have developed a multi-scale mathematical model of the female rainbow trout hypothalamus-pituitary-ovary-liver axis to use as a tool to help understand the functioning of the system and for extrapolation of laboratory findings of stressor impacts on specific components of the axis. The model describes the essential endocrine components of the female rainbow trout reproductive axis. The model also describes the stage specific growth of maturing oocytes within the ovary and permits the presence of sub-populations of oocytes at different stages of development. Model formulation and parametrization was largely based on previously published in vivo and in vitro data in rainbow trout and new data on the synthesis of gonadotropins in the pituitary. Model predictions were validated against several previously published data sets for annual changes in gonadotropins and estradiol in rainbow trout. Estimates of select model parameters can be obtained from in vitro assays using either quantitative (direct estimation of rate constants) or qualitative (relative change from control values) approaches. This is an important aspect of mathematical models as in vitro, cell-based assays are expected to provide the bulk of experimental data for future risk assessments and will require quantitative physiological models to extrapolate across biological scales.


Introduction
Fishes, like other vertebrates, have an endocrine signaling network operating between the brain and the gonads controlling reproduction. This network spans the hypothalamus, the pituitary, the gonads and, in fishes and other egg-laying vertebrates, the liver, and is referred to in this paper as the hypothalamus-pituitary-ovary-liver (HPOL) axis (Fig 1 and Table 1) [1]. Communication between these components takes place via the blood as well as through neurosecretory fibers from the hypothalamus to the pituitary [2]. The HPOL axis begins in the hypothalamus where environmental cues trigger the production of gonadotropin-releasing hormone (GnRH) [3]. GnRH is then sent to the pituitary where it stimulates the synthesis of two gonadotropins (GTHs): follicle-stimulating hormone (FSH) and luteinizing hormone (LH) [4]. Both GTHs regulate production of sex steroids by ovarian follicles, where the follicles are defined as the oocyte and the surrounding cell layers [5,6]. FSH stimulates the somatic cell layers of the follicle to produce estradiol-17β (E2) [7], which stimulate the liver to produce the egg-yolk protein vitellogenin (VTG). VTG is subsequently secreted into the bloodstream and incorporated into the developing oocyte, the accumulation of which is the primary cause for oocyte growth. LH is involved in the later stages of oocyte development and stimulates the somatic cell layers of the follicles to produce a maturation-inducing steroid, which in rainbow trout (Onchorynchus mykiss) and many other fishes is 17α, 20β-dihydroxy-4-pregnen-3-one (DHP). DHP is the primary hormone involved in the final stages of oocyte maturation and ovulation [8]. Beyond these functions, the above five hormones interact with each other through positive and negative feedback loops, creating an intricate network of pathways that leads to the synchronization of processes culminating in the timely production of mature eggs.
The overall complexity of the fish reproductive axis is on a level similar to mammalian systems; however, development of mathematical descriptions for fish has only occurred relatively recently. One challenge with modeling the fish reproductive axis is that spawning behavior generally occurs in two forms: group synchronous spawning (a single large clutch of oocytes develop synchronously for one spawning event; e.g. salmon, trout) or asynchronous spawning (several small clutches of eggs are spawned at different times during a reproductive season; e.g. zebrafish [Danio rerio], fathead minnow [Pimephales promelas], medaka [Oryzias latipes]). Although the same set of reproductive hormones appear to be involved in both types of fishes, the roles of some hormones, such as the GTHs, are not as well defined in asynchronous compared to synchronous spawning fishes. This and other factors associated with oocyte development have made conceptualization of models distinct between the two spawning types. In group synchronous spawning fishes an early pituitary-gonad-liver model for the Atlantic croaker (Micropogonias undulates) was described [9,10] and included descriptions of production for a single GTH (combined FSH\LH), E2 formation and subsequent VTG production. The model was then used to predict the effects of chemical exposures (a PCB mixture and cadmium) and environmental stress (prolonged hypoxia) on the pituitary-gonad-liver axis [9,10]. A model for the coho salmon (O. kisutch) hypothalamus-pituitary-gonad axis was described [11], which The HPOL signaling network in rainbow trout as formulated in our model. Arrows and symbols on graph follow CellDesigner vs. 4.4 notation (www.celldesigner.org). GnRH is secreted from the hypothalamus into the pituitary stimulating the production of mFSH and mLH, which then leads to formation of FSH and LH, respectively. FSH, which is being continuously secreted from the pituitary, travels to the ovaries to stimulate production of E2. E2 then travels to the liver to bind with E2 receptors (R; translated from mR) to form ER. ER then stimulates the production of mVTG, which produces VTG L . Secreted VTG then travels from the liver to the ovaries via the plasma (VTG P ) where it is absorbed by follicles in stages 3 through 6 (the proportion of follicles in these stages are denoted by S j , j = 3, 4, 5, and 6) during vitellogenesis, the rate of which is affected by FSH P , to promote oocyte growth (O Avg ). Oocyte growth then progresses the oocytes through the stages using a Weibull distribution created from O Avg together with O Var . In the later stages LH P stimulates the oocytes to produce DHP. Finally, oocytes undergo final maturation (S FOM ) and combined with DHP, determine when the fish ovulates. Description of symbols are found in Table 1. added a simple description of oocyte maturation, synthesis of DHP and separate descriptions for FSH and LH. A more detailed description of vitellogenesis throughout the reproductive cycle was described for rainbow trout [12] and included parameters for the expression of mRNA coding for the estrogen receptor and VTG. A specialized model for the unnatural stimulation of vitellogenesis in male trout was presented [13], which focused on transcriptional imprinting of VTG expression. In asynchronous spawning fish, greatest attention has been placed on the fathead minnow and a focus towards its importance as an ecotoxicological test species. A series of pituitary-gonad-liver models for this species have been described [14,15,16,17]. These models included empirical descriptions of a single GTH variable and subsequent stimulation of testosterone and E2 synthesis in the gonads. Model formulation was integrated with toxicokinetic descriptions of E2 or ethynylestradiol (EE2; synthetic E2 agonist) exposure to male minnows [14], EE2 and 17β-trenbolone (synthetic testosterone agonist) to female minnows [16] or fadrozole (E2 synthesis inhibitor) to female minnows [15,17]. These models were used to characterize the reproductive effects of exposure to sex steroid agonists or antagonists to fathead minnows. Other more specialized models for the fathead minnow include the model by [18] for E2 synthesis, which characterizes the multiple steps involved in the conversion of cholesterol to E2, and the model by [19], which models the process of asynchronous oocyte development by tracking multiple clutches of oocytes through a maturation cycle and spawning. The emphasis in modeling E2 and VTG production in fathead minnows reflects the critical role these variables play in controlling fecundity and successful spawning [20].
The growing interest in modeling the fish reproductive axis is due in part to awareness that accelerating climate change, habitat destruction and discharging of contaminants into waterways have potential effects on natural and cultured fish populations. New approaches to understanding the impacts of these stressors are being developed that require more biologically accurate and detailed mathematical models [21,22]. The value of these models for mammalian systems has been well-demonstrated and resides in their potential to aid in understanding biological control processes associated with reproduction and test hypotheses that may be difficult or impractical to do in vivo [23]. The objective of the present study is to present a comprehensive model for the HPOL axis for female rainbow trout that can be used to assess impacts of physical and chemical stressors on reproduction. We combine elements of past models [11,12] with new features that better characterize processes associated with oocyte growth and maturation.

Ethics statement
All experiments were performed according to the guidelines established by the Institutional Animal Care and Use Committee (IACUC) of PNNL.

Background
The central focus of the model is on the progression of events that lead to oocyte growth through several stages, followed by final oocyte maturation (FOM) and then ovulation, whose timing is taken here to be the time point at which DHP first exceeds a threshold and when almost all of the oocytes are in FOM. We do not include actual spawning in the model because the release of eggs is different in the wild (where it is triggered by interactions with male fish) and in aquaculture (where the eggs are usually "stripped" by hand). In the model, concentrations of FSH, LH and VTG are given for separate tissue compartments. To differentiate between compartments subscripts are used: blood plasma is denoted with a P, the pituitary is denoted with a Pit, the ovary with a O and the liver is denoted with an L. Other tissues not specified are denoted with an N. The description of the equations for a hormone's plasma levels occurs in the compartment associated with its synthesis and secretion. A clearance-volume approach is used to describe the pharmacokinetic behavior of VTG and circulating hormones, with a plasma referenced volume of distribution and a clearance parameter that represents loss due to excretion or sequestration (e.g. VTG accumulation in oocytes). The volume of distribution of a hormone is denoted by V i and plasma clearance of a hormone is denoted by Cl i , where i is replaced with the hormone name, e.g. FSH, LH, E2. In addition to the previously mentioned parameters, a consistent notation is used for parameters performing the same biological function, e.g. synthesis and degradation, to help maintain model uniformity.

Functions
The model uses Hill and Heaviside (a.k.a. the Unit Step function) functions to describe nonlinear hormone interactions. The Hill function is a sigmoidal function with values ranging from 0 to 1 and is used to model stimulatory (H + ) and inhibitory (H -) effects that have a threshold, H þ ðx; T; nÞ ¼ x n x n þ T n and H À ðx; T; nÞ ¼ In each equation x is the concentration of the hormone, protein, or mRNA, T is the threshold value at which x exerts half of its maximal feedback and n determines the range of concentrations x must fall in to switch between no effect and maximal effect. Values of x not close to T will fall outside of the interval T±ε n , where ε n is a positive constant depending on n and is a mathematical representation of the maximal distance a value can be from T and still be considered close. For values of x outside this interval the Hill function remains relatively close to zero, no effect, or relatively close to one, maximal effect. As n becomes larger ε n becomes smaller and the change between no effect and the maximal effect becomes more sudden; letting n go to infinity results in an instantaneous switch between no effect, zero, and the maximal effect, one, when concentration levels of x reach T. This on/off switch can be modeled using the Heaviside function, by evaluating it at x-T when letting n go to infinity in H + or T-x when letting n go to infinity in H -.
To incorporate time delays in interactions between model variables, transit compartments are used (reviewed in [24]). Transit compartments have been shown to be well suited to account for delays associated with signal transduction processes [25], which for example, occurs in the present model with FSH stimulation of E2 synthesis. This approach allows for a less complex model but still retains the capacity to incorporate more details of the signal transduction process, which can directly replace transit compartments. An individual transit compartment is restricted to short time delays, so for longer time delays that indicate multiple distinct processes are occurring, more than one transit compartment is used. Furthermore, the number of distinct processes can increase or decrease based on level of detail describing the biological process. This allows us to not only find the optimal number of transit compartments to use to represent the delay of hormone A's influence on hormone B by D A,B hours, denoted by m, but it is also reasonable to assume the delay produced by an individual transit compartment is equal to the average delay of the compartments, D A,B /m. A superscript of TC m is used for a system of m transit compartments; for example, if the growth of B starting at time 0 is modeled by then incorporating the time delay using a system of m transit compartments will change the formula to where the transit compartments are given by The use of transit compartments, Hill functions and the Heaviside function help describe delays and nonlinear interactions between hormones and key reproductive processes.

Hypothalamus-Pituitary
The hypothalamic region of the mid-brain is well established as the signal initiation site that stimulates the pituitary gland at the base of the brain to release GTHs that act on the gonads to cause their sexual development. The endocrine control of the HPOL pathway begins in neurosecretory fibers from the hypothalamus that project directly into the adenohypophysis region of the pituitary gland [26]. These fibers release GnRH in the vicinity of the gonadotrophs, the pituitary cells that produce the GTHs, increasing the basal synthesis rate of the beta subunit mRNA for FSH (mFSH) and LH (mLH). The GTHs are heterodimeric proteins consisting of a common α and GTH specific β subunit.
The regulation of GnRH is partly due to kisspeptins (KISS), which are the links from environmental cues to the HPOL axis [2]. The complexity of the GnRH-KISS system is still being characterized in fish and therefore we have chosen to not include it in this version of the model. We have assumed that FSH beta subunit mRNA levels primarily reflect changes in GnRH released into the pituitary during the reproductive cycle. Thus, an empirical equation for GnRH production was derived based on previously measured mFSH levels in the trout pituitary. The equations for mFSH and mLH are based on the equations found in [11] and begin with a basal synthesis rate, k s,mFSH and k s,mLH , and a degradation rate, k d,mFSH and k d, mLH . The positive feedback of GnRH is represented through a stimulatory factor, α mFSH,GnRH and α mLH,GnRH . E2 is known to increase the expression of mLH and pituitary content of LH [27,28,29]. This positive feedback effect of E2 on mLH is represented though a stimulatory factor, α mLH,E2 , and delayed through the use of transit compartments by D E2,mLH hours to account for signal transduction and other biological processes.
After FSH is synthesized from mFSH it is immediately released into the blood. The tight coupling of mFSH levels and blood FSH levels ( [30], and data presented in this study) allows the model to have a constant rate of synthesis and release of FSH from the pituitary, k s,FSH , adjusted by the weight of the pituitary, w Pit . As in [11], the clearance-volume approach then leads to the following equation: Plasma levels of FSH are detectible throughout the reproductive cycle, implying a continuous release from the pituitary; in contrast, plasma levels of LH remain largely undetectable until the end of the reproductive cycle. Though plasma levels of LH remain undetectable, pituitary levels continue to increase [30] implying a block on the release of LH from the pituitary into the plasma. The release of LH from the pituitary is known to be controlled by a dopamine (D2) inhibition, which is strengthened by E2 through positive feedback on D2 receptors (D2r) in the pituitary [31,32]. Declining levels of E2 remove the block on LH release [32] only after levels surpass a threshold, T E2,LH , [33]. Besides E2, rising DHP levels also appear to block LH release. In some fishes, DHP causes an increase of D2 receptors in the pituitary [31], and in rainbow trout, DHP significantly decreased the release of LH in-vitro [34]. To model the release of LH from the pituitary we have chosen a release function that assumes the following: for each ng/ml of LH in the pituitary there exists a minimum amount of D2 receptors needed to maintain the block, the amount of LH released depends on the pituitary levels of LH that do not have the necessary amount of D2r to be blocked, the quantity of D2 receptors is proportional to levels of E2 and DHP and for the block to be removed, E2 levels must be higher than a threshold T E2,LH . To better understand how these assumptions lead to the final release function found in Eq (12) an initial release function will be created from the first assumption, Eq (9), and built upon from the remaining assumptions, Eq (10) through Eq (12). The first assumption can be interpreted as the release of LH occurs when the amount of LH in the pituitary exceeds the minimum amount of D2-receptors per ng/ml of LH in the pituitary, ξ where ξ>0, e.g.
Applying the second and third assumption to the release function gives where . Finally, the fourth assumption gives the final release function of With the exception of the release function the equations modeling the pituitary and plasma levels of LH are similar to those found in [11]. In the pituitary LH is synthesized from mLH at a constant rate of k s,LH . The stored LH is then subject to degradation at a rate of k d,LH before it is released at a rate of k r,LH . The LH released into the plasma is adjusted by the weight of the pituitary, w Pit , and follows the clearance-volume approach.
After release into the plasma, the GTHs travel to the ovaries to stimulate production of the sex steroids.

Ovary
The process of oocyte growth can be initially separated according to the capacity to absorb vitellogenin: pre-vitellogenic, vitellogenic and post-vitellogenic. These stages can be further subdivided according to specific morphological and biochemical changes occurring during growth [6]. Our model divides the reproductive cycle into seven stages; stages 1 and 2 are pre-vitellogenic, stages 3 through 6 are vitellogenic and the last stage is FOM. The proportion of oocytes in stage j, j = 1,2,. . .,6,FOM, at time t will be denoted by S j (t). The model in [11] characterizes the stage of an oocyte based on length of time in the reproductive process and did not contain feedback from other variables. In the present model the stage of an oocyte is based on its size, which is defined by its diameter, and interacts with other variables. This is a convenient metric that is commonly reported in experimental studies and has been used to determine oocyte staging in trout (see Appendix S1 for ways to translate diameter to other measurements of size). The distribution of diameters at time t is given by a Weibull distribution with time dependent parameters l(t) and k(t). The distribution's parameters are chosen such that the mean and variance of the Weibull distribution at time t is equal to the mean and variance of oocyte sizes given by Eqs (16) and (18), respectively. Hence the stages are given by Pðs jÀ1 X t s j Þ j ¼ 2; . . . ; 6 where s j is the boundary size between stages j-1 and j and the random variable X t is the size of an oocyte at time t with X t * Weibull(l(t),k(t)). When describing oocyte growth it is important to recognize that approximately 84% of the increase in oocyte diameter occurs during vitellogenesis [5]. Therefore, the model assumes that oocyte growth rate is dependent on the amount of VTG being sequestered, Seq(t), and the conversion of diameter increase for every unit of VTG sequestered,k V;O Avg . The growth rate in the non-vitellogenic stages is described by a constant growth rate,k NV;O Avg , as only a small portion of growth occurs during these stages.
The rate at which the oocytes absorb VTG is determined in part by circulating levels of FSH, oocyte stage, and VTG levels [35,36,37]. The increase in basal sequester rate, Cl VTG,seq , due to circulating FSH is represented with a Hill function with a half threshold of T Seq,FSH . Since oocytes can only sequester VTG when they exceed a certain size, [37], and have not entered FOM, the sequester period will only be from stage 3 until stage 6.
The growth rate of individual oocytes can vary substantially, presumably due to differential uptake of VTG but other, unknown factors may also be involved. At the start of the reproductive cycle, all primary oocytes (stage 1) appear to be uniform in size. As oocytes mature towards the mid-vitellogenic stage, size variance increases before becoming more uniform as final maturation nears [38]. To model this, each stage will have a maximal variance assigned to it,a O Var ;S j , where the variances increase as the oocytes near the mid-vitellogenic stage, Stage 4, and decrease as the oocytes approach FOM.
During stages 2 through FOM oocytes produce E2 and DHP when FSH and LH, respectively, are introduced to the ovaries. Once in the ovaries FSH and LH stimulate the thecal and granulosa cells, the layers of somatic cells surrounding the oocytes, to produce E2 and DHP, respectively [1]. When FSH or LH binds to their respective receptors, the ensuing stimulation of steroid synthesis is not instantaneous. To account for this delay, transit compartments are introduced in place of FSH and LH in Eqs (19) and (20), respectively, with time delays D FSH,E2 and D LH,DHP , respectively. The sensitivity of the oocytes' responses to LH and FSH is stage specific and the amount of steroid secreted into the plasma can be thought of as a clearance from the follicle. Hence a follicle in stage j will produce Cl E2;S j (Cl DHP;S j ) ml of E2 (DHP) for each ng/ ml of FSH (LH) that stimulates. Therefore, the total amount of E2 produced during Stage j at time t is Cl E2;S j Á n oocyte Á S j ðtÞ Á FSH TC n P Â Ã , where n oocyte is the average number of oocytes per kg of fish. Similarly, the total amount of DHP produced during Stage j at time t is In addition to FSH stimulated E2 production, we have a basal production rate for each follicle, k E2 . We do not include a basal production rate for DHP since the follicles ability to produce DHP is only a short time period.
Once produced, E2 and DHP are then released into the blood and stimulate the production of VTG and promote ovulation, respectively.

Liver
The liver plays an essential role in the reproductive process because it synthesizes the egg yolk precursor protein VTG, which is necessary for oocyte growth and development [39]. Our model for VTG synthesis is largely based on the model described by [12]. VTG synthesis is initiated when E2 binds to a nuclear estrogen receptor, R, at a binding rate of k on,ER , which is produced by transcription of the estrogen receptor mRNA, mR. The resulting E2-receptor complex, ER, stimulates the production of mR in a positive auto-regulatory feedback loop and stimulates the production of the mRNA for VTG, mVTG. Based on the experimental data of [40] we assume that activation of transcription happens on a one-to-one ratio, meaning a given dose of E2 yields similar increases in the E2-receptor complex, ER, and mVTG levels. ER then dissociates at a rate of k off,ER increasing the number of unbound E2-receptors, R.
Once synthesized, each molecule of mVTG can be translated into VTG multiple times; this ability is represented by the amplification factor γ. Furthermore, N mVTG represents a scaling factor for mVTG concentration. The VTG in the liver, denoted by VTG L , is secreted into the plasma at a rate of k r,VTG , the amount of which is denoted by VTG P . The amount of VTG released into the plasma is measured in g/kg of liver weight, w L . VTG is then exchanged between the plasma and compartments other than the liver and ovaries, VTG N , taken up by the ovaries through the oocytes, Eq (17), or cleared from the body, Cl VTG . The VTG in compartments other than the ovary are assumed to be in reversible equilibrium with blood with a transfer clearance of Cl VTG,trans .
Parameter identification, model calibration and sensitivity analysis Most parameter values were estimated directly from previously published experimental studies, prior modeling efforts [11,12] or modified slightly from literature values. These parameters and their sources are identified in Table 2. Other parameters were estimated using observed data collected from female rainbow trout during a second reproductive cycle. The latter data was previously published in [41] (oocyte diameter; plasma E2 and VTG; liver ER and VTG mRNA) or measured for this study (plasma FSH, LH, DHP and pituitary mRNA for beta subunits of FSH and LH). The data for FSH and LH was analyzed by radioimmunoassay as described by [42]. DHP was determined by enzyme immunoassay using reagents purchased from Cayman Chemical (Ann Arbor, MI). The β subunit mRNA for FSH was measured by q-RT-PCR as described in [43]. The β subunit mRNA for LH was also measured by q-RT-PCR using a similar method except the forward and reverse primers were from the LH-β gene (Gen-Bank accession number AB050836). These primer sequences were: forward-AAAC-GATCCGCCTACCTGACT and reverse-AGCCACAGGGTAGGTGACATG. From these data, input curves were created using MATLAB's cubic spline program and used to decouple the system of ODEs into smaller subsystems with fewer unknown parameters. Using the remaining input curves and the solution of the subsystem of ODEs, the error was then optimized around a least squares cost function with respect to the unknown parameters. After various decouplings, the solution for the original system of ODEs was compared to the input equations using a least squares cost function. The final cost function was then optimized by making small perturbations in all parameters not obtained explicitly through the literature. The parameters gained from this optimization process are also in Table 2. After the parameters were obtained, a global sensitivity analysis was performed using the PAWN method, outlined in [44], where the cumulative distribution functions created by a one-dimensional output distribution are used to compute density-based sensitivity indices. The model predictions of mFSH, mLH, FSH P , LH Pit , LH P , O Avg , E2, DHP, mR, R, ER, mVTG, VTG L , VTG P and VTG N using the original parameters found in Table 2 and a set of altered parameters were compared using the Euclidean norm to measure the relative changes at every hour over a 350 day time period. Summing the relative changes of the individual functions provides a one-dimensional representation of the HPOL model that normalizes the effects the altered parameters have on the individual functions. The set of altered parameters was created using a Sobol' quasi-random sequence to obtain a uniform set of parameters that vary between ±20% of the original values listed in Table 2. For the HPOL model varying the parameters between ±20% of the original values provided a large diverse set of model predictions while lowering the proportion of  similar results. These results are displayed in Table 2 with larger numbers indicating a more sensitive parameter that exerts greater influence. Many of the more sensitive parameters are those associated with the synthesis and pharmacokinetics of E2 and VTG and growth of oocytes (Table 2; highest sensitivity values). Within the pituitary, the most sensitive parameters are associated with FSH synthesis. The greater sensitivity of these parameters reflects the model's emphasis on oocyte growth and staging. This corresponds with the biology in that oocyte growth over the entire cycle is driven by FSH through its effect on E2, which stimulates VTG production.

Results
The lifespan of rainbow trout in their natural habitat is quite variable and may last from three to seven years [54]. Most populations of rainbow trout mature in their second or third year Parameter value was obtained directly from source. 2 Guided initial estimates of the parameter value. 3 Parameter value was obtained from prior modeling efforts. 4 Specified data was used to create an input curve and an initial estimate was obtained through optimization methods solving a minimum number of differential equations. 5 Initial estimate was obtained by calculating time delays in the relative extrema from specified data. 6 Final parameter value was refined through optimization methods using experimental data from [41]  and may spawn in successive years up to five times, typically in late fall or early spring [55,56,57]. A complete reproductive cycle in female rainbow trout normally requires 345-365 days, which is reproduced in the model. We've chosen to focus on changes in model variables during a second reproductive cycle because a large portion of the experimental data that was used to guide model development was collected during this time period. Model simulations were performed using ode23s in MATLAB (version R2014a) and the parameter values are listed in Table 2. A summary of model predictions for tissue specific and circulating variables are presented in Figs 2-4. Included in these figures are experimentally measured values from a cohort of females sampled during a second reproductive cycle. Model simulations were compared to this data set as it represents, to our knowledge, the most complete set of measurements of model variables in a single cohort of female rainbow trout. Although portions of this data set were used to guide estimation for some parameters (notably synthesis rates of mFSH and mLH), many others were determined from separate in vivo and in vitro studies (Table 2) and the model was not optimized or fitted to this data set. The model is also capable of continuous predictions, which we demonstrate for circulating variables (FSH, LH, E2, DHP and VTG) over three successive spawning cycles (Fig 2).  Table 2. Observed data (mean ± SD, n = 3) measured from a cohort of second time spawning female trout is shown as dark grey circles. The observed data for FSH, LH and DHP was measured as part of this study or from [45]. Three weeks after spawning average oocyte growth was reset to O Avg (0). Circulating DHP is shown on a nonlinear scale to emphasize low levels (less than 1 ng/ml) of DHP before FOM. Figs 2 and 3 suggest the model adequately describes the changes that occur during the reproductive cycle. Model predicted levels of circulating FSH reflect the pattern of mFSH synthesis in the pituitary throughout the cycle and is in close agreement with observed data (Figs 2A and 3A). After reaching the ovaries, FSH stimulates the follicles to produce E2 (Fig 2C), the rate of which is stage dependent with increasing rates until FOM (Table 2). Hence, the relative extrema of circulating E2 levels reflect the relative extrema in circulating levels of FSH after an approximately 36 day delay. This activity of FSH is known to be mediated by a membrane bound receptor in the thecal and granulosa cells of the follicle [1]. We have not included GTH receptor dynamics in this version of the model and instead describe the stimulation of E2 synthesis with a single rate parameter.

Inspection of model predicted and observed values in
Although not modeled explicitly, we assume that a portion of the stage specific changes in the rate of E2 synthesis is associated with changes in expression of the FSH receptor. When circulating levels of E2 increase beyond a threshold, subsequent declines trigger the release of LH.  Table 2. Observed data (dark grey circles; mean ±TG mRn = 3) was measured from the same individuals used for Fig 2. doi:10.1371/journal.pcbi.1004874.g003 The model predicts a premature rise in circulating LH that occurs approximately 80 days before the main surge in LH that triggers the final changes leading to ovulation (Fig 2B). This early or premature rise in LH is not reflected in the observed data although it may have been easily missed due to the rapid increase and decline of LH in blood. One possible role for an early LH increase is to stimulate DHP synthesis, which has been suggested to promote meiosis in early oogonia that will become oocytes during the next cycle [58]. The model over predicts mVTG synthesis in the liver (Fig 3D) although prediction of circulating VTG closely matches observed data (Fig 2E).
As the oocytes progress into stage 6 and FOM, rising DHP levels eventually pass a threshold and fulfill one of two requirements for ovulation. The other requirement for ovulation depends on the proportion of oocytes that reach FOM, which in the model depends on the growth and variance of the oocytes (Fig 4A and 4B respectively). Initially, the distribution of oocytes is positively skewed with little variation in size. As vitellogenesis begins the growth rate of individual oocytes varies depending on a number of factors causing a large variation in sizes; the peak of which occurs during mid-vitellogenesis, stage 4. When the oocytes enter the post-vitellogenic phases and near FOM, the growth rate becomes inversely related to size causing a reduction in variation and the distribution to be negatively skewed. This gradual increase and then steep decline in oocyte variance (as seen in Fig 4D) causes the early and mid vitellogenic stages to be the largest growth stages (Fig 4C). Subsequent rapid growth of oocytes during stages 5 and 6, days 225 through 335, is caused by the steep increase in VTG production ( Fig 4A). When the oocytes enter FOM they lose the capacity to sequester VTG causing a transient spike in circulating VTG levels (Fig 2E).
Because parameterization of the model relied on multiple in vitro and in vivo sources of experimental data, it is reasonable to expect the model to be relatively robust in describing reproductive performance beyond the data set used for initial validation (Figs 2-4). To demonstrate this, we present two examples using experimental data from three additional sources. These studies monitored several model variables (FSH, LH, VTG, E2) in cohorts of secondtime spawning rainbow trout [53,59,60]. Comparison of HPOL model predictions with observed data from these studies is shown on Figs 5 and 6. In both examples, subtle differences in observed and predicted FSH levels lead to differences in either E2 (Fig 5; [59]) or VTG (Fig 6; [60]). However, the model is still able to reasonably describe oocyte growth patterns from both studies (Figs 5 and 6). If the model's GnRH function is optimized for the FSH profiles for each study, an improvement in the prediction of E2 and VTG is observed (green lines on Figs 5 and 6). By adjusting the oocyte sizes that define stages 3-6 to account for approximately 10%, 9%, 31% and 34% of oocyte growth (black lines on Figs 5 and 6), model predictions of oocyte growth more closely match the observed values. These adjustments made to the oocyte stages follow a similar logic to the derivation of the parameter dividing stage 6 from FOM, s 6 , in Table 2 and uses the parameter values listed in Table 2 to approximate the percentage of growth each of the vitellogenic stages should account for. These examples highlight the ability of the model to predict reproductive variables in rainbow trout but also expose some limitations of broadly using the model with rainbow trout as currently parameterized. There is considerable plasticity in spawning behavior among rainbow trout, which may spawn at various times throughout the year [57]. Thus, it should be expected that some natural variability in GnRH and FSH synthesis will occur in different fish populations and this variation will introduce changes in model variables associated with E2 and VTG synthesis. However, the generally good agreement of model simulations with previously published data supports validation of the model and the parameters gained from optimization in Table 2 and when necessary, should require only minor adjustments to describe reproduction in different populations of rainbow trout.

Discussion
Reproduction in fishes and other vertebrates represents the timely coordination of many endocrine factors that culminate in the production of mature, viable gametes. The model presented here accurately describes the essential processes of GTH production, oocyte growth, vitellogenesis and final maturation in a group synchronous annual spawning fish. The model is an improvement over previous efforts in salmonids, which focused more narrowly on either GTH and steroid production [11] or on vitellogenesis [12]. Here we have combined and expanded these models into a single version that can continuously describe reproduction in trout over successive cycles. The present model provides a description of GnRH production and its stimulation of GTH synthesis throughout the reproductive cycle. This is in contrast to [11] where GnRH was modeled at a constant level of 50 arbitrary units between days 75 and 250 of the reproductive cycle. This simplified approach proved inadequate to capture the seasonal changes observed for pituitary expression of mFSH and secretion of FSH into plasma. Accurate prediction of mFSH expression is particularly important due to the tight coupling with  [59]. Model simulations using parameter estimates listed in Table 2 and the same GnRH function used in Figs 2-4 are shown as dashed blue lines. Simulations using a customized GnRH function based on measured FSH levels in [59] (see also Appendix S5) are shown as a green dashed line. Simulations with the customized GnRH function and altered vitellogenic staging parameters (s 3 = 1.19, s 4 = 1.74, s 5 = 3.51 and s 6 = 5.48) are shown as a black solid line. A regression curve was used to map the GSI data in [59] to the average oocyte diameter (See Appendix S1). Another specific improvement in the present model is the more extensive description of oocyte staging and oocyte growth. The model explicitly describes oocyte size, which is used to define growth and to differentiate the developmental stages of maturation. Including distinct  [53,60]. Model simulations using parameter estimates listed in Table 2 and the same GnRH function used in Figs 2-4 are shown as dashed blue lines. Simulations using a customized GnRH function based on measured FSH levels in [53,60] (see also Appendix S5) are shown as a green dashed line. Simulations with the customized GnRH function and altered vitellogenic staging parameters (s 3 = 1.21, s 4 = 1.76, s 5 = 3.48 and s 6 = 5.39) are shown as a black solid line. A regression curve was used to map the GSI data in [53] to the average oocyte diameter (See Appendix S1). stages of development is advantageous because it is a traditional measurement reported by fish biologists and many processes such as VTG accumulation are stage dependent. Also, within the maturing trout ovary there is evidence for sub-populations of oocytes with overlapping developmental stages and associated stage dependent growth rates [38]. Defining the stages provides the ability to include attributes that are stage specific and allow for overlap. The majority of oocytes will fall into one stage category at a time, allowing for segregation of stages. However, segregated staging may cause discontinuities in the ODE system and would not satisfy the criteria for uniqueness and existence of solutions for ODEs, the Picard-Lindelöf theorem [61]. Although smooth approximations could be used to define the staging to avoid discontinuities, the solution would be irregular (See Appendix S4). An added attraction of overlapping staging is that it narrows the distinction made between synchronous and asynchronous spawning fish and offers the potential for future studies to explore control processes that contribute to the evolution of reproductive strategies.
A challenging aspect of HPOL model development is how to approach the timely release of LH from the pituitary into the bloodstream near the end of the reproductive cycle. In constructing the model, we evaluated several equations of varying complexity for describing the release of LH (see Appendix S2). The release primarily depends on declining E2 levels that typically occur during the later stages of oocyte development. Hence, LH release could simply be restricted to occur when E2 declines during a narrowly defined time period such as during late vitellogenesis and FOM: This approach is similar to the release function used in [11]. However, there are reports of elevated blood levels of LH before the final decline of E2 [41,59,62] and other evidence suggesting the release of LH is not confined to the final stages of reproduction [33]. Therefore, a more biologically realistic approach would be to remove restrictions on stage specificity and allow LH release to occur when E2 levels are declining faster than a threshold (T E2 ) . Another consideration is whether secretion of LH from the pituitary can simply be treated as a bolus release (e.g. near instantaneous), which has the effect of causing a large spike in predicted blood levels of LH. Although adequate to describe the LH plasma profile, this also tends to overestimate peak levels of LH because the actual release probably occurs as a more gradual process. It also ignores the role that dopamine plays in maintaining the block on LH release. This led us to consider the release of LH as the spillover of unblocked LH in the pituitary, Eq (12). Alternatively, a more physiologically explicit description of dopamine signaling and levels of D2r could be added to the model. However, we have chosen not to add this to the model at this time because of the added complexity and lack of experimental data in trout to validate model predictions.
Another important aspect of model development was to incorporate time delays associated with biological processes not specifically accounted for in the model. An example of this is the receptor-mediated effects of FSH on the follicle, particularly with respect to steroidogenesis, which occur over an extended time period. Rather than attempt to include additional parameters to characterize FSH activity, we chose to use transit compartments to represent time delays that are due to receptor signaling, gene expression and associated changes in biochemical processes. An alternative to transit compartments would be the use of delayed differential equations (DDEs). The advantages of DDEs reside with the potential to describe more complex behavior and improve accuracy with less parameters [24]. However, the use of DDEs requires knowledge of a variable's behavior before the initial time point, information which may not be available or difficult to ascertain during initial ovarian recrudescence. Transit compartments are a simple ODE approximation of DDEs which only require knowledge of present conditions and are a well-established tool in pharmacokinetic-pharmacodynamic modeling to characterize time delays [25]. Furthermore, transit compartments can be readily exchanged with submodels to provide a more detailed approach to specific biological processes. For example, in the present model E2 is synthesized and added to the bloodstream after FSH enters the ovaries with a time delay of D E2,FSH hours. A more detailed model of the synthesis of E2 could be used such as the type described by [18] and inserted in place of the transit compartments between FSH and E2. Conversely, the present model could be simplified by inserting transit compartments between E2 and plasma VTG to approximate plasma level of VTG (see Appendix S3 for results): This latter approach might be used when the model is being applied to other fish species that lack detailed experimental data on vitellogenesis. In either case, the use of transit compartments provides flexibility for increasing or decreasing complexity as needed.

Application and future work
An important value of quantitative biological models is the ability to perform simulations that vary model parameters in ways that are associated with various physiological conditions, real or hypothetical. In this manner, the rainbow trout HPOL model can be a tool to explore the effects of diverse physical and chemical stressors on reproduction. As an example, we consider trenbolone exposure to female rainbow trout. Trenbolone is a synthetic testosterone analogue that is used to promote growth in cattle. Concerns over the possible release of trenbolone to aquatic environments have prompted several studies on its effects to fish. In fathead minnows, trenbolone exposures (27 ng/L and higher) for up to three weeks while fish are actively spawning, caused a rapid and significant decline in circulating E2 and VTG levels and eventual cessation of spawning [51,63,64]. Subsequent in vitro experiments with fathead minnow ovary tissue suggest the decrease was due in part to decreased expression of enzymes involved in E2 synthesis [65]. In short-term exposures (one-four weeks) to female rainbow trout, trenbolone exposure stimulated FSH release in pre-vitellogenic fish, but also caused decreases in E2 levels in both pre-and mid-vitellogenic trout [51]. Other experiments with trout indicated plasma clearance of E2 was increased by 50% [51]. We were interested in estimating whether these effects in rainbow trout alone or in various combinations would be sufficient to prevent reproduction if trout were assumed to be continuously exposed throughout the reproductive cycle. We adjusted the following parameters based on the prior experimental evidence: FSH release (k s,FSH ; increased 4x), E2 synthesis (SE2O4-5; decreased 46%), E2 plasma clearance (Cl E2 ; increased 41%) and VTG expression (k s,mVTG ; decreased 40%). A summary of model simulations for E2 and oocyte diameter is shown in Fig 7. Increasing the release of FSH alone or in combination with up to 46% decreases in E2 synthesis did not affect overall oocyte growth (Fig  7). If a decrease in VTG expression or increased E2 clearance occurs simultaneously with decreased E2 synthesis, then oocyte growth is slowed to an extent that ovulation would be unlikely to occur (Fig 7). These simulations suggest trenbolone effects on E2 synthesis by itself are not sufficient to cause reproductive failure in trout and additional effects occurring simultaneously elsewhere in the system are needed. Additional experimental research is needed to verify these observations but the exercise illustrates the potential for the HPOL model to incorporate in vivo and in vitro results for predictions of reproductive performance. This is a valuable attribute because long-term or complete reproductive cycle exposures are costly to perform and require many animals.
An emerging concept in environmental health research is the adverse outcome pathway (AOP) [66]. An AOP seeks to extrapolate effects of stressors measured at the cellular / sub-cellular level of biological organization to whole organism and population-level effects. It is anticipated that future research efforts will increasingly rely on in vitro methods to generate health effects data. Thus, mathematical models such as the trout HPOL model presented here will be useful for translating in vitro based observations into whole fish effects. As we consider future efforts to develop the next-generation HPOL model, emphasis will be placed on linking model outputs more closely to the types of input data used for population models. This would include expanding the description of oocyte development to provide more detail on egg quality and lifetime fecundity of rainbow trout over multiple spawning cycles.
Supporting Information S1 Fig. Predicted plasma levels of LH using different release functions over an entire reproductive cycle, (A), and from day 250 until the end of the reproductive cycle, (B). The black dashed line is using the final release function, R LH , given by Eq (14). The red line is using the release function R LH1 given by Eq (28). The blue line is using the release function R LH2 given by Equation (S5). The green line is using the release function R LH3 given by Equation (S7). Plasma LH data from second time spawning female rainbow trout, also used in Fig 2,  Model predicted profiles of plasma levels of (A) E2 and (B) average oocyte growth for rainbow trout exposed to trenbolone. Previously published studies suggest trenbolone exposure may decrease E2 synthesis to varying degrees depending on exposure levels, and also increase the clearance of E2 from plasma and increase the secretion of FSH in trout. The solid black line is the model predictions for fish not exposed to trenbolone and uses the parameter values found in Table 2. The dashed lines depict model predictions of known effects of trout that are exposed to trenbolone and use parameter values found in Table 2 except for Cl E2,Sj = 0.54*Cl E2,Sj for j = 4 and 5, k s,FSH = 4*k s,FSH during Stages 1 and 2. The blue dashed line also uses Cl E2 = 1.41*Cl E2 and the green dashed line uses k s,mVTG = 0.6*k s,mVTG . The model predicts that ovulation is unlikely to occur when an increase in E2 clearance and/or a decrease in mVTG synthesis is added to the effects because the delay in ovulation is more than 50 days.