# An optimization formulation for characterization of pulsatile cortisol secretion

^{1}Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, USA^{2}Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA^{3}Department of Anesthesia, Critical Care and Pain Medicine, Massachusetts General Hospital, Boston, MA, USA^{4}Laboratory for Information and Decision Systems, Massachusetts Institute of Technology, Cambridge, MA, USA^{5}Engineering Systems Division, Massachusetts Institute of Technology, Cambridge, MA, USA^{6}Institute for Data, Systems, and Society, Massachusetts Institute of Technology, Cambridge, MA, USA^{7}Institute for Medical Engineering and Science, Massachusetts Institute of Technology, Cambridge, MA, USA^{8}Department of Anesthesia, Harvard Medical School, Boston, MA, USA

Cortisol is released to relay information to cells to regulate metabolism and reaction to stress and inflammation. In particular, cortisol is released in the form of pulsatile signals. This low-energy method of signaling seems to be more efficient than continuous signaling. We hypothesize that there is a controller in the anterior pituitary that leads to pulsatile release of cortisol, and propose a mathematical formulation for such controller, which leads to impulse control as opposed to continuous control. We postulate that this controller is minimizing the number of secretory events that result in cortisol secretion, which is a way of minimizing the energy required for cortisol secretion; this controller maintains the blood cortisol levels within a specific circadian range while complying with the first order dynamics underlying cortisol secretion. We use an ℓ_{0}-norm cost function for this controller, and solve a reweighed ℓ_{1}-norm minimization algorithm for obtaining the solution to this optimization problem. We use four examples to illustrate the performance of this approach: (i) a toy problem that achieves impulse control, (ii) two examples that achieve physiologically plausible pulsatile cortisol release, (iii) an example where the number of pulses is not within the physiologically plausible range for healthy subjects while the cortisol levels are within the desired range. This novel approach results in impulse control where the impulses and the obtained blood cortisol levels have a circadian rhythm and an ultradian rhythm that are in agreement with the known physiology of cortisol secretion. The proposed formulation is a first step in developing intermittent controllers for curing cortisol deficiency. This type of bio-inspired pulse controllers can be employed for designing non-continuous controllers in brain-machine interface design for neuroscience applications.

## 1. Introduction

Many hormones that have been well-investigated appear to be released in pulses (Stavreva et al., 2009); for example, cortisol, gonadal steroids, and insulin are released in a pulsatile manner (Veldhuis, 2008). Pulsatility is a physiological way of increasing hormone concentrations rapidly and sending distinct signaling information to target cells (Veldhuis, 2008). Ultradian pulsatile hormone secretion allows for encoding information via both amplitude and frequency modulation and is a way of frequency encoding (Lightman and Conway-Campbell, 2010; Walker et al., 2010b). Pulsatile signaling permits target receptor recovery, rapid changes in hormone concentration, and greater control, and is also more efficient than continuous signaling (Walker et al., 2010b). The mechanism underlying the generation of hormone pulses and why this method of signaling is chosen by the body over continuous signaling is not known. Since the transcriptional program prompted by hormone pulses is considerably different from constant hormone treatment (Stavreva et al., 2009), it is crucial to understand the physiology underlying pulsatile hormone release. Hormone pulsatility underlies multiple physiological processes. For example, (i) cortisol oscillations have crucial effects on target cell gene expression and glucocorticoids receptor function (McMaster et al., 2011; Walker et al., 2012). (ii) Some psychiatric and metabolic diseases are associated with changes in cortisol pulsatility (Walker et al., 2010a). (iii) When the same amount of corticosterone is administered by constant infusion rather than a pulsatile infusion, it results in a noticeably reduced ACTH response to stress (Lightman and Conway-Campbell, 2010). In this study, we investigate pulsatile release of cortisol and propose a novel mathematical formulation that characterizes pulsatile cortisol secretion.

Cortisol is released from the adrenal glands in pulses in response to pulsatile release of ACTH. CRH induces the release of ACTH. In return, cortisol has a negative feedback effect on ACTH and CRH release at the pituitary and hypothalamic levels. The timing and amplitudes of cortisol pulses vary throughout the day where the amplitude variations are due to the circadian rhythm underlying cortisol release with periods of 12 and 24 h (Faghih et al., 2011), and the variations in the timing of cortisol pulses result from the ultradian rhythm underlying cortisol release. Between 15 and 22 secretory pulses of cortisol are expected over 24 h (Veldhuis et al., 1989; Brown et al., 2001).

Based on the interactions in the HPA axis, it was hypothesized that pulsatile release of CRH from the hypothalamus results in pulsatile release of cortisol. Walker et al. suggest that a sub-hypothalamic pituitary-adrenal system results in the pulsatile ultradian pattern underlying cortisol release (Walker et al., 2012). This is because inducing constant CRH levels results in a pulsatile cortisol profile (Walker et al., 2012) while constant ACTH levels do not result in pulsatile cortisol secretion (Spiga et al., 2011). Spiga et al. suppressed the activity of the HPA axis by oral methylprednisolone and infused both constant amounts and pulses of ACTH to test the hypothesis that pulsatile ACTH release is necessary for pulsatile cortisol secretion (Spiga et al., 2011). While pulsatile ACTH resulted in pulsatile cortisol secretion, constant infusion of same amounts of ACTH did not activate cortisol secretion (Spiga et al., 2011). Moreover, studies on sheep in which the hypothalamus has been disconnected from the pituitary suggest that pulsatile input from hypothalamic secretagogues (e.g., CRH or vasopressin) is not necessary for the ultradian rhythm in cortisol secretion or for pulsatile cortisol secretion and pulsatile cortisol secretion is still maintained (Walker et al., 2010a). Hence, pulsatile cortisol release is controlled by the dynamics in the anterior pituitary. Since pulsatile cortisol release seems to be more efficient than continuous signaling, it might be the case that the anterior pituitary is solving an optimal control problem.

We postulate that there is a controller in the anterior pituitary that controls the pulsatile secretion of cortisol and the ultradian rhythm of the pulses via the negative feedback effect of cortisol on the anterior pituitary. Hence, by considering the known physiology of the HPA axis, we shall formulate an optimization problem that achieves impulse control. In optimal control theory, impulse control is a special case of bang-bang control, in which an action leads to instantaneous changes in the states of the system (Sethi and Thompson, 2000). Impulse control occurs when there is not an upper bound on the control variable and an infinite control is exerted on a state variable in order to cause a finite jump (Sethi and Thompson, 2000). Minimizing an ℓ_{0}-norm cost function can achieve impulse control and we use a reweighed ℓ_{1}-norm formulation as a relaxation to the ℓ_{0}-norm to solve the proposed optimization formulation. Moreover, we consider the first-order dynamics underlying cortisol synthesis and the circadian amplitude constraints on the cortisol levels when formulating the optimization problem.

## 2. Methods

We propose a physiologically plausible optimization problem for cortisol secretion by making the following assumptions: (1) Cortisol levels can be described by first-order kinetics for cortisol synthesis in the adrenal glands, cortisol infusion to the blood, and cortisol clearance by the liver described in Brown et al. (2001), Faghih (2010), and Faghih et al. (2011, 2014). (2) There is a time-varying cortisol demand [*h*(*t*)] that should be satisfied throughout the day, which is a function of the circadian rhythm. (3) There is a time-varying upper bound on the cortisol level [*q*(*t*)] that is a function of the upper bound on the cortisol level that the body can produce or a holding cost so that the cortisol level would not be much above the demand. (4) Control that results in cortisol secretion [*u*(*t*)] is non-negative. (5) The body is minimizing the number of resources (control) throughout the day. Hence, we postulate that there is a controller in the anterior pituitary that controls cortisol secretion via the following optimization formulation:

where *x*_{1} is the cortisol concentration in the adrenal glands and *x*_{2} is the blood cortisol concentration. λ and γ, respectively, represent the infusion rate of cortisol from the adrenal glands into the blood and the clearance rate of cortisol by the liver.

Considering the known physiology of *de novo* cortisol synthesis (i.e., no cortisol is stored in the adrenal glands) (Brown et al., 2001), we assume that the initial condition of the cortisol level in the adrenal glands is zero [*x*_{1}(0) = 0] (Brown et al., 2001). Assuming that the input and the states are constant over 1-min intervals, and *y*_{0} is the initial condition of the blood cortisol concentration, blood cortisol levels at every minute over *N* min can be represented in discrete form by ${\text{y}}{=}{{\left[}\begin{array}{cccc}{{y}}_{{1}}& {{y}}_{{2}}& {\cdots}& {{y}}_{{N}}\end{array}{\right]}}^{{\prime}}$ where *y*_{k} is the blood cortisol level at time *k* and ${\text{y}}$ can be represented as:

where ${F}{=}{{\left[}\begin{array}{cccc}{{f}}_{{1}}& {{f}}_{{2}}& {\cdots}& {{f}}_{{N}}\end{array}{\right]}}^{{\prime}}{,}{{f}}_{{k}}{=}{{e}}^{{-}{\gamma}{k}}{,}{G}{=}{{\left[}\begin{array}{cccc}{{g}}_{{1}}& {{g}}_{{2}}& {\cdots}& {{g}}_{{N}}\end{array}{\right]}}^{{\prime}}$, ${{g}}_{{k}}{=}{{\left[}\begin{array}{ccc}\frac{{\lambda}}{{\lambda}{-}{\gamma}}{(}{{e}}^{{-}{\gamma}{k}}{-}{{e}}^{{-}{\lambda}{k}}{)}& {\mathrm{...}}& \frac{{\lambda}}{{\lambda}{-}{\gamma}}{(}{{e}}^{{-}{\gamma}}{-}{{e}}^{{\lambda}}{)}\underset{{N}{-}{k}}{\underbrace{{0}{}{\cdots}{}{0}}}\end{array}{\right]}}^{{\prime}}$, and ${\text{u}}$ represents the control over *N* min. Then by letting ${\text{h}}{=}{{\left[}\begin{array}{cccc}{{h}}_{{1}}& {{h}}_{{2}}& {\cdots}& {{h}}_{{N}}\end{array}{\right]}}^{{\prime}}$ where *h*_{k} is the cortisol demand at an integer minute *k* and ${\text{q}}{=}{{\left[}\begin{array}{cccc}{{q}}_{{1}}& {{q}}_{{2}}& {\cdots}& {{q}}_{{N}}\end{array}{\right]}}^{{\prime}}$ where *q*_{k} is the upper bound at the integer minute *k*. Hence, we solve the discrete analog of the formulation in Equation (1):

ℓ_{0} problems are generally NP-hard, and instead an ℓ_{1}-norm relaxation of such problems can be solved. In solving ℓ_{1}-norm problems, there is a dependence on the amplitude of the coefficients over which the ℓ_{1}-norm is minimized, and there is more penalty on larger coefficients than on smaller ones. However, it is possible to strategically construct a reweighted ℓ_{1}-norm such that non-zero coefficients are penalized in a way that the cost further resembles the ℓ_{0}-norm. By putting large weights on small entries, the solution concentrates on entries with small weights, non-zero entries are discouraged in the recovered signal, and a cost function that is more similar to an ℓ_{0}-norm cost function can be solved (Candes et al., 2008). To find such weights for ℓ_{1}-norm cost function, Candes et al. (2008) have proposed an iterative algorithm for enhancing the sparsity using reweighted ℓ_{1} minimization, which solves $\underset{{\text{u}}}{{\mathrm{min}}}{\Vert}{\text{u}}{{\Vert}}_{{0}}$. This algorithm is based on Fazel's “log-det heuristic” algorithm for minimizing the number of non-zero entries of a vector (Fazel, 2002) and the convergence of this log-det heuristic algorithm has been studied in Lobo et al. (2007). Hence, we use the algorithm by Candes et al. (2008) such that the constraints in the optimization problem in Equation (3) are satisfied:

1. Initialize the diagonal matrix *W*^{(0)} with entries *w*^{(0)}_{i} = 1, *i* = 1, …, *n* on the diagonal and zeros elsewhere

2. Solve ${{\text{u}}}^{{(}{\ell}{)}}{=}{a}{r}{g}\underset{{\text{u}}}{{\text{}}{\mathrm{min}}}{\Vert}{{W}}^{{(}{\ell}{)}}{\text{u}}{{\Vert}}_{{1}}$ subject to the constraints in Equation (3)

3. Update the weights ${{w}}_{{i}}^{{(}{\ell}{+}{1}{)}}{=}\frac{{1}}{{\left|}{{\text{u}}}_{{i}}{}^{{(}{\ell}{)}}{\right|}{+}{\u03f5}}$, *i* = 1, …, *n*

4. Iterate till ℓ reaches a certain number of iterations. Otherwise, increment ℓ and go to step 2.

The idea is that by solving ${{\text{u}}}^{{(}{\ell}{+}{1}{)}}{=}{\text{}}{\text{arg}}\underset{{\text{u}}}{{\mathrm{min}}}{\displaystyle {{\sum}}_{{i}{=}{1}}^{{n}}\frac{{\left|}{{\text{u}}}_{{i}}{\right|}}{{\left|}{{\text{u}}}_{{i}}{}^{{(}{\ell}{)}}{\right|}{+}{\u03f5}}}$ iteratively, the algorithm attempts to solve for a local minimum of a concave penalty function that is more similar to the ℓ_{0}-norm (Candes et al., 2008). ϵ is used to ensure that weights on the recovered zero entries will not be set to ∞ at the next step, which would prevent us from obtaining estimates at the next step. ϵ should be slightly larger than the expected non-zero amplitudes of the signal that is to be recovered, and a value of at least 0.001 is recommended (Candes et al., 2008). This algorithm does not always find the global minimum and as ϵ → 0, the likelihood of stagnating at an undesirable local minimum increases (Candes et al., 2008). For ϵ values closer to zero, the iterative reweighted ℓ_{1}-norm algorithm stagnates at an undesirable local minimum (Candes et al., 2008).

We study the optimization problem in Equation (1) via four examples. We first investigate the case that the optimization formulation in Equation (1) is selecting the control such that the state (i.e., the blood cortisol concentration) is bounded between constant lower and upper bounds to illustrate the idea that the formulation in Equation (1) can achieve impulse control. Then, we investigate cases in which the upper and lower bounds have harmonic profiles with a circadian rhythm. Using the iterative algorithm for enhancing the sparsity by reweighted ℓ_{1} minimization (Candes et al., 2008), we solve the optimization problem in Equation (1) over a time period τ and update the solution after a time period $\frac{{\tau}}{{2}}$ and repeat this process for a 24-h period. λ, γ, ϵ, τ, and lower and upper bounds are given in Tables 1–3. Since empirically the algorithm converges in 10 iterations for the formulation in this study, we use ℓ = 10 when running the algorithm. Numerical analysis was performed in MATLAB R2011b and using CVX (Grant and Boyd, 2008, 2014).

## 3. Results

To illustrate that the proposed approach results in impulse control, we use constant lower and upper bounds and show that the proposed method achieves impulse control and a state that has a pulsatile profile. This example is not physiological and is used to help the reader better understand the type of results this type of approach generates. Then, we show an example that corresponds to a healthy subject and leads to impulse control. The secretory events and cortisol levels are in agreement with physiologically plausible profiles in healthy human data, and the obtained solution is optimal. Moreover, we illustrate another example that corresponds to a healthy subject and achieves impulse control. In this example, while the secretory events and cortisol levels are physiologically plausible, the obtained solution is optimal over the first 20 h and suboptimal for the last 4 h. This example shows that the performance of the algorithm used for solving the proposed optimization formulation depends on the choice of ϵ and can stagnate at a local minimum. Finally, we provide an example that illustrates a case in which the number of pulses is not within a physiologically plausible range (i.e., an abnormality) while impulse control is achieved.

### 3.1. Example 1

Assuming that the upper and lower bounds are constant, the optimal solution is achieved when the initial condition starts at the upper bound; then, when the state decays to the lower bound, an impulse causes a jump in the state which brings it back to the upper bound, and then again the state decays to the lower bound and the same jump to the upper bound again occurs, and the same process keeps repeating. Figure 1 shows that solving the optimization problem (Equation 1) for constant upper and lower bounds using the parameters given for Example 1 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively, results in impulse control. There are 12 constant impulses obtained over a 24-h period, which occur periodically. This example is just a simple toy problem illustrating that the optimization formulation in Equation (1) can achieve impulse control and pulsatile cortisol release using a low energy input. This example does not have any physiological implications for cortisol secretion as it does not include upper and lower bounds that have a circadian rhythm observed in cortisol levels.

**Figure 1. Cortisol levels and control obtained using Example 1**. **(A)** The top panel displays the optimal cortisol profile (black curve), constant upper bound (red curve), and constant lower bound (blue curve). **(B)** The bottom panel displays the optimal control. The optimization problem obtained 12 impulses over 24 h as the optimal control (the timing of the control was discretized into 1440 points; the obtained control takes 12 non-zero values, i.e., impulses, while it is zero everywhere else). The optimization problem was solved using the parameters given in Example 1 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively.

### 3.2. Example 2

In healthy humans, cortisol levels have regular periodic time-varying patterns that consist of episodic release of secretory events with varying timings and amplitudes in a regular diurnal pattern. Figure 2 shows that solving the optimization problem (Equation 1) for two-harmonic bounds with a circadian rhythm, using the parameters given for Example 2 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively, the obtained control is impulse control. Figure 2 also displays that adding a zero mean Gaussian measurement error with a standard deviation of σ = 0.45 to each simulated data point and recording the cortisol levels every 10 min (which is comparable to measurement noise and sampling interval of cortisol data in human subjects, Faghih et al., 2014), the obtained cortisol profile resembles cortisol human data provided in Faghih et al. (2014). There are 16 impulses over a 24-h period with time-varying circadian amplitudes and ultradian timings; the obtained control is within the physiologically plausible range of 15 –22 pulses (Veldhuis et al., 1989; Brown et al., 2001). The impulses are more frequent during the day and have higher amplitudes during the day than in night time. Obtained cortisol levels are low at night. Then, around 6 AM, cortisol levels increase, reaching higher values between 10 AM and 12 PM, followed by a gradual decrease throughout the day reaching low values at night. The obtained control and state are optimal; the state starts at the upper bound and decays to the lower bound at which point an impulse causes a jump in the system that results in increasing the state, and the state reaches the upper bound. Then, the state decays again to the time-varying lower bound and this process repeats. This example illustrates that the optimization formulation in Equation (1) can achieve impulse control and pulsatile cortisol release, using a low energy input, and generate secretory events and cortisol levels that have physiologically plausible profiles similar to those observed in healthy human data.

**Figure 2. Cortisol levels and control obtained using Example 2**. **(A)** The top panel displays the optimal cortisol profile obtained by adding a zero mean Gaussian measurement error with a standard deviation of σ = 0.45 to each simulated data point; the cortisol levels are recorded every 10 min. **(B)** The middle panel displays the optimal cortisol profile (black curve), two-harmonic upper bound (red curve), and two-harmonic lower bound (blue curve); the cortisol levels are recorded every minute. **(C)** The bottom panel displays the optimal control. The optimization problem obtained 16 impulses over 24 h as the optimal control (the timing of the control was discretized into 1440 points; the obtained control takes 16 non-zero values, i.e., impulses, while it is zero everywhere else). The optimization problem was solved using the parameters given in Example 2 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively.

### 3.3. Example 3

In this example, we consider different lower and upper bounds compared to Example 2 while keeping λ and γ to values used in Example 2. Figure 3 shows that solving the optimization problem (Equation 1) for two-harmonic bounds with a circadian rhythm, using the parameters given for Example 3 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively, the obtained control is impulse control. Figure 3 also displays that adding a zero mean Gaussian measurement error with a standard deviation of σ = 0.45 to each simulated data point and recording the cortisol levels every 10 min (which is comparable to measurement noise and sampling interval of cortisol data in human subjects, Faghih et al., 2014), the obtained cortisol profile resembles cortisol human data provided in Faghih et al. (2014). Sixteen impulses are obtained over 24 h which is within the physiological range of 15–22; these impulses have time-varying circadian amplitudes and ultradian timings. The impulses have higher amplitudes and are more frequent between 4 AM and 12 PM. The obtained cortisol levels are low at night. Then, the cortisol levels increase, reaching higher values between 7 AM and 11 AM, followed by a gradual decrease throughout the day, reaching low values at night. This example illustrates that the optimization formulation in Equation (1) can achieve impulse control and pulsatile cortisol release using a low energy input, and generates secretory events and cortisol levels that have physiologically plausible profiles similar to those observed in healthy human data. The control and state obtained in the first 20 h are optimal; however, the control and the state obtained for the last 4 h are suboptimal as the algorithm used for solving the optimization problem (Equation 1) can stagnate at a local minimum depending on the choice of ϵ. However, still a low energy control is recovered that keeps the cortisol levels within the desired bounds.

**Figure 3. Cortisol levels and control obtained using Example 3**. **(A)** The top panel displays the cortisol profile obtained by adding a zero mean Gaussian measurement error with a standard deviation of σ = 0.45 to each simulated data point; the cortisol levels are recorded every 10 min. **(B)** The middle panel displays the obtained cortisol profile (black curve), two-harmonic upper bound (red curve), and two-harmonic lower bound (blue curve). **(C)** The bottom panel displays the obtained control. The optimization problem obtained 16 impulses over 24 h as the control (the timing of the control was discretized into 1440 points; the obtained control takes 16 non-zero values, i.e., impulses, while it is zero everywhere else). The optimization problem was solved using the parameters given in Example 3 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively.

### 3.4. Example 4

In this example, we keep the lower and upper bounds the same as the values we used in Example 3 while using values for λ and γ that result in higher infusion of cortisol and lower clearance of cortisol compared to Example 3. Figure 4 shows that solving the optimization problem (Equation 1) using the parameters given for Example 4 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively, the obtained control is impulse control. Figure 4 also displays that adding a zero mean Gaussian measurement error with a standard deviation of σ = 0.45 to each simulated data point and recording the cortisol levels every 10 min (which is comparable to measurement noise and sampling interval of cortisol data in human subjects Faghih et al., 2014), the obtained cortisol profile resembles cortisol human data provided in Faghih et al. (2014). Twelve impulses are obtained over 24 h where the impulses have lower amplitudes and are less frequent compared to the impulses obtained in Example 3. The obtained impulses still have time-varying circadian amplitudes and ultradian timings. The number of pulses has decreased compared to Example 3 which was expected as cortisol is cleared faster in this example. While the number of these pulses are not within the physiological range reported for healthy subjects, the obtained cortisol levels are still within the desired range. Cortisol levels are low at night, then increase, reaching higher values between 6 AM and 10 AM, followed by a gradual decrease throughout the day, reaching low values at night. The peak values of cortisol levels change and on average in this example the cortisol levels have lower values, and this might illustrate a case of cortisol deficiency. Also, in this example, the optimization formulation in Equation (1) results in impulse control and pulsatile cortisol release using a low energy input. The control and state obtained in the first 19 h are optimal; however, the control and the state obtained for the last 5 h are suboptimal as the algorithm used for solving the optimization problem (Equation 1) can stagnate at a local minimum depending on the choice of ϵ. However, still a low energy control is recovered that keeps the cortisol levels within the desired bounds.

**Figure 4. Cortisol levels and control obtained using Example 4**. **(A)** The top panel displays the cortisol profile obtained by adding a zero mean Gaussian measurement error with a standard deviation of σ = 0.45 to each simulated data point; the cortisol levels are recorded every 10 min. **(B)** The middle panel displays the obtained cortisol profile (black curve), two-harmonic upper bound (red curve), and two-harmonic lower bound (blue curve). **(C)** The bottom panel displays the obtained control. The optimization problem obtained 12 impulses over 24 h as the control (the timing of the control was discretized into 1440 points; the obtained control takes 12 non-zero values, i.e., impulses, while it is zero everywhere else). The optimization problem was solved using the parameters given in Example 4 in Table 1 and the upper and lower bounds provided in Tables 2, 3, respectively.

## 4. Discussion

It is well-known that cortisol is released in pulses, and based on our results it appears that this method of relaying information might be an optimal approach as opposed to continuous signaling. In this work, we formalized this concept by proposing an optimization formulation for a physiologically plausible controller in the anterior pituitary that achieves impulse control as the optimal solution. In the proposed formulation, we assumed that there is a time-varying upper bound on the cortisol levels in the blood. Also, we assumed that the cortisol levels in the blood should be above a time-varying circadian threshold to achieve normal regulation of the HPA axis. We assumed that the lower bound and upper bound on the cortisol levels are two-harmonic functions with periods of 12 and 24 h that are controlled by the circadian rhythm. However, the upper bound and the lower bound for cortisol secretion could have multiple harmonics, and this assumption is only considering the most significant periods in cortisol release. Moreover, we considered the first-order dynamics underlying cortisol secretion. We have shown that the proposed optimization formulation yields impulse control as its optimal solution. The number, timing, and amplitude of the recovered secretory events in the proposed optimization problem are physiologically plausible. Moreover, the obtained cortisol profile is in agreement with the circadian rhythm observed in healthy human data. As pointed out, the iterative algorithm for enhancing the sparsity by reweighted ℓ_{1} minimization (Candes et al., 2008) does not always find the global minimum and might stagnate at an undesirable local minimum; we employed this algorithm to solve examples of optimization problems formulated in Equation (1) to show that the formulation in Equation (1) achieves impulse control as observed in cortisol levels. However, the optimization problem in Equation (1) can be solved using other methods as well, and for arbitrary choices of ϵ and τ the algorithm for enhancing the sparsity by reweighted ℓ_{1} minimization (Candes et al., 2008) might stagnate at a local minima and not achieve the optimal solution (please see Example 3).

To validate this mathematical characterization using experiments, one can start by recovering the parameters for a rat model and obtain lower and upper bounds on cortisol levels in a healthy rat. Next, make the adrenal glands of the rat malfunctional so that the rat becomes Addisonian and does not secrete cortisol. Then, using a pulse controller, obtain a cortisol profile that stays within the lower and upper bounds found when the rat was healthy.

While we proposed a simple optimization formulation that can achieve impulse control, it is possible to obtain impulse control using more complex formulations by either assuming that the system is a switched system with different rates or assuming that the nature of the system is impulsive and there is no continuous control. We assumed that the infusion and clearance rates are constant; however, the system can be a switched system with different infusion and clearance rates. Abrupt changes in the infusion and clearance rates could also result in impulse control. For example, if the infusion rate of cortisol from the adrenal glands starts from a constant level at wake and decreases abruptly to a new constant level, a very large level of cortisol should be produced in a short time so that the desired cortisol level can still be achieved. There could be multiple abrupt changes in the infusion rate throughout the day, and there might be an infusion rate reset to a high level at the beginning of sleep. Another example that could possibly result in impulse control is when the clearance starts at a constant level, and increases abruptly to a new constant level; then, a very large level of cortisol should be produced in a short time so that the desired cortisol level can still be achieved. There could be multiple such abrupt changes in the clearance rate throughout the day, and the clearance rate might be reset to a low level at the beginning of sleep. Another scenario could be that both the infusion and the clearance rates could be starting from a constant level and change abruptly to different levels periodically. In that case, the overall effect is that cortisol gets cleared faster or cortisol gets infused to the blood more slowly, and at such moments a very large cortisol level should be released for a short period of time to maintain the desired cortisol level. Such situations could possibly achieve impulse control as long as there is not an upper bound on the control variable; a mathematical example of a model with a time-varying rate that achieves impulse control is given in Sethi and Thompson (2000), and the *maximum principle* is used to find the optimality conditions for this problem. Moreover, it is possible that pulsatile inputs arise from the nature of the system, and the hormone system might be designed such that the input to the system can only be impulsive where the timing of the impulses are functions of the states and are not activated until a reseting condition is satisfied. A mathematical example of such a model is given in Wang and Balakrishnan (2008) where the cost function minimizes the energy in the input and the state, and calculus of variations is used to find the optimality conditions. Also, another possibility is that the body is solving a weighted ℓ_{1} cost function where different costs are associated with the control at different times of the day (e.g., the weights obtained at convergence when using the reweighted algorithm).

In this study, for modeling cortisol secretion, we proposed a physiologically plausible optimization formulation for a controller in the anterior pituitary. A similar approach can be used to study other endocrine hormones that are released in pulses. For example, the proposed optimization formulation can be tailored to include the constraints underlying thyroid hormone secretion or gonadal hormone secretion or growth hormone secretion in order to study the pulsatile release of those hormones. The transcriptional program stimulated by hormone pulses is very different from constant hormone treatment and some disorders are associated with hormone pulsatility. Hence, understanding the underlying nature of the pulsatile release of these hormones via mathematical formalization can be beneficial to understanding the pathological neuroendocrine states and treating some hormonal disorders.

In addition to contributing to the scientific advances in understanding cortisol regulation in daily rhythms, we provide a better understanding of the biological mechanism mathematically, which can potentially be used to come up with a similar approach to devise pulsatile control interventions instead of continuous controllers for treating cortisol disorders. Traditional control-theoretic methods do not normally consider the intermittent control that is observed in pulsatile control of cortisol release. Instead of developing a controller that tracks the desired cortisol levels, we have proposed a formulation for a controller that maintains the cortisol levels within certain upper and lower bounds. Our study formalizes, mathematically, the pulsatile controller underlying cortisol secretion, and through a simulation study we show that our formulation can control the cortisol levels to remain within the desired bounds while having the circadian and the ultradian rhythms underlying cortisol secretion. Hence, while our approach uses control-theoretic concepts to understand a biological process, the proposed formulation is a first step in developing intermittent controllers for curing cortisol deficiency. While the methods proposed in our study are not externally applied to control a biological process, with slight modifications based on the pathological condition of interest, the proposed intermittent controller can be used to control some of the pathological problems related to cortisol. This can be done by including the first-order kinetics of the medicine that will be injected to the patient to control cortisol levels, and then using compressed sensing algorithms to recover the secretory release of cortisol in the patient. In this case there will be two sets of pulses that control cortisol levels: (i) external pulses that are injected to the patient (ii) pulses that are secreted as a part of the natural control system underlying cortisol secretion. Similarly, such bio-inspired controllers can be used for controlling other hormones (e.g., growth hormone, thyroid hormone, or gonadal hormones).

Since this type of controllers can be adapted to be applied for curing different pathological conditions related to endocrine hormones, the idea behind modeling such controllers opens new research directions. For example, a patient who suffers from Addisons disease takes cortisone once or twice a day for their cortisol deficiency (which does not seem optimal), while an impulse controller can be used to control the cortisol levels optimally. The future directions of this research include designing an impulse controller such that the optimality of the controller is guaranteed. Moreover, in brain-machine interface design, in which brain implants control epilepsy or Parkinson's disease, it is possible to design pulse controllers instead of continuous controllers to improve the battery life of the brain implant and reduce the number of surgeries required for changing the battery of the implanted controller. With the new advances and ongoing research in brain-machine interface design for psychiatric disorders, this type of pulse controller can potentially be used to control post-traumatic stress disorder, major depression, and addiction. For example, in post-traumatic stress disorder or major depression, in theory, one could potentially measure skin conductance response which results from discrete emotional shocks experienced by the patient, and ideally stimulate ventromedial prefrontal cortex using impulse control to reverse the effect of the emotional shocks in the patient. In conclusion, inspired by the pulse controller proposed in this research, it is potentially possible to design a class of pulse controllers for applications that naturally arise in neuroscience.

## Author Contributions

RTF, MAD, and ENB designed the optimal control formulation. RTF performed research and wrote the paper.

## Funding

RTF's work was supported in part by the NSF Graduate Fellowship. For this work, ENB is supported in part by NIH DP1 OD003646, 1-R01-GM104948-03 and NSF 0836720, and MAD is supported in part by EFRI-0735956. The funders had no role in study design and analysis, decision to publish, or preparation of the manuscript.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

RTF would like to thank Dr. Gail K. Adler for her feedback on endocrinology. This work was presented in RTF's Ph.D. thesis (Faghih, 2014).

## References

Brown, E. N., Meehan, P. M., and Dempster, A. P. (2001). A stochastic differential equation model of diurnal cortisol patterns. *Am. J. Physiol. Endocrinol. Metab.* 280, E450–E461. Available online at: http://ajpendo.physiology.org/content/280/3/E450.long

Candes, E. J., Wakin, M. B., and Boyd, S. P. (2008). Enhancing sparsity by reweighted ℓ_{1} minimization. *J. Fourier Anal. Appl.* 14, 877–905. doi: 10.1007/s00041-008-9045-x

Faghih, R. T. (2010). *The Fitzhugh-Nagumo Model Dynamics with an Application to the Hypothalamic Pituitary Adrenal Axis*. S.M. thesis, Massachusetts Institute of Technology, Cambridge.

Faghih, R. T. (2014). *System Identification of Cortisol Secretion: Characterizing Pulsatile Dynamics*. Ph.D. thesis, Massachusetts Institute of Technology, Cambridge.

Faghih, R. T., Dahleh, M. A., Adler, G. K., Klerman, E. B., and Brown, E. N. (2014). Deconvolution of serum cortisol levels by using compressed sensing. *PLoS ONE* 9:e85204. doi: 10.1371/journal.pone.0085204

Faghih, R. T., Savla, K., Dahleh, M. A., and Brown, E. N. (2011). “A feedback control model for cortisol secretion,” in *2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society* (Boston, MA), 716–719.

Fazel, M. (2002). *Matrix Rank Minimization with Applications*. Ph.D. thesis, Stanford, CA: Stanford University.

Grant, M., and Boyd, S. (2008). “Graph implementations for non-smooth convex programs,” in *Recent Advances in Learning and Control, Lecture Notes in Control and Information Sciences*, eds V. Blondel, S. Boyd, and H. Kimura (Springer-Verlag Limited), 95–110. Available online at: http://stanford.edu/~boyd/graph_dcp.html

Grant, M., and Boyd, S. (2014). *CVX: Matlab Software for Disciplined Convex Programming, Version 2.1.* Available online at: http://cvxr.com/cvx

Lightman, S. L., and Conway-Campbell, B. L. (2010). The crucial role of pulsatile activity of the HPA axis for continuous dynamic equilibration. *Nat. Rev. Neurosci.* 11, 710–718. doi: 10.1038/nrn2914

Lobo, M. S., Fazel, M., and Boyd, S. (2007). Portfolio optimization with linear and fixed transaction costs. *Ann. Oper. Res.* 152, 341–365. doi: 10.1007/s10479-006-0145-1

McMaster, A., Jangani, M., Sommer, P., Han, N., Brass, A., Beesley, S., et al. (2011). Ultradian cortisol pulsatility encodes a distinct, biologically important signal. *PLoS ONE* 6:e15766. doi: 10.1371/journal.pone.0015766

Sethi, S. P., and Thompson, G. L. (2000). *Optimal Control Theory: Applications to Management Science and Economics*. Boston, MA: Kluwer Academic Publishers.

Spiga, F., Waite, E. J., Liu, Y., Kershaw, Y. M., Aguilera, G., and Lightman, S. L. (2011). ACTH-dependent ultradian rhythm of corticosterone secretion. *Endocrinology* 152, 1448–1457. doi: 10.1210/en.2010-1209

Stavreva, D. A., Wiench, M., John, S., Conway-Campbell, B. L., McKenna, M. A., Pooley, J. R., et al. (2009). Ultradian hormone stimulation induces glucocorticoid receptor-mediated pulses of gene transcription. *Nat. Cell Biol.* 11, 1093–1102. doi: 10.1038/ncb1922

Veldhuis, J. D. (2008). “Pulsatile hormone secretion: mechanisms, significance and evaluation,” in *Ultradian Rhythms from Molecules to Mind*, eds D. Lloyd and E. L.Rossi (Springer), 229–248.

Veldhuis, J. D., Iranmanesh, A., Lizarralde, G., and Johnson, M. L. (1989). Amplitude modulation of a burstlike mode of cortisol secretion subserves the circadian glucocorticoid rhythm. *Am. J. Physiol. Endocrinol. Metab.* 257, E6–E14.

Walker, J. J., Spiga, F., Waite, E., Zhao, Z., Kershaw, Y., Terry, J. R., et al. (2012). The origin of glucocorticoid hormone oscillations. *PLoS Biol.* 10:e1001341. doi: 10.1371/journal.pbio.1001341

Walker, J. J., Terry, J. R., and Lightman, S. L. (2010a). Origin of ultradian pulsatility in the hypothalamic–pituitary–adrenal axis. *Proc. R. Soc. Lond. B Biol. Sci.* 277, 1627–1633. doi: 10.1098/rspb.2009.2148

Walker, J. J., Terry, J. R., Tsaneva-Atanasova, K., Armstrong, S. P., McArdle, C. A., and Lightman, S. L. (2010b). Encoding and decoding mechanisms of pulsatile hormone secretion. *J. Neuroendocrinol.* 22, 1226–1238. doi: 10.1111/j.1365-2826.2010.02087.x

Keywords: pulsatile control, cortisol secretion, endocrine control, mathematical modeling, circadian rhythm

Citation: Faghih RT, Dahleh MA and Brown EN (2015) An optimization formulation for characterization of pulsatile cortisol secretion. *Front. Neurosci*. 9:228. doi: 10.3389/fnins.2015.00228

Received: 29 March 2015; Accepted: 11 May 2015;

Published: 11 August 2015.

Edited by:

Jason Ritt, Boston University, USAReviewed by:

Alireza Mousavi, Brunel University, UKJose Luis Contreras-Vidal, University of Houston, USA

Copyright © 2015 Faghih, Dahleh and Brown. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Rose T. Faghih, Neuroscience Statistics Research Laboratory, Massachusetts Institute of Technology, 77 Massachusetts Avenue, 46-6057 Cambridge, MA, USA, rfaghih@mit.edu