A Computer Simulation Model for Automated Quantification of Luteinizing Hormone Secretion

: Problem statement: An automated software package is developed for quantifying Luteinizing Hormone (LH) release from pituitary gonadotropes in response to short pulses of gonadotropin-releasing hormone (GnRH). Approach: These computer programs were designed to accommodate the release mechanism that consists of three stages, namely; (1) Production of different concentrations of effectors (such as receptor-hormone binds, G protein interactions, inositol 1,4,5-trisphosphate IP 3 ), (2) Release of Ca 2+ from the endoplasmic reticulum ER and finally (3) Release of LH through pumping of Ca 2+ in and out of the cytosol. Results: The programs were coded in C++ computer language and implemented on Pentium PC. The designed model was intended to help in explaining the characteristics of LH release in response to various duration and different concentrations of GnRH pulses in the presence and absence of external Ca 2+ . Conclusion: Features of the short-term responses were simulated and fairly understood using receptor binding, dimerization, interaction of the dimerized receptor with a G protein, production of effectors that open voltage sensitive channels in cell membrane and finally open the Ca 2+ channels in the ER and the Ca 2+ - dependent release of LH. Also, prediction of the varying concentrations of Ca 2+ effects and comparison with experimental data was made possible.


INTRODUCTION
The tremendous increase in the computer speed, memory size, variety of visualization software techniques and hardware components, besides the growing decrease in hardware cost have made it the most attractive tool for various applications [1,2] . Utilization of computer abilities in the medical field has no limit, ranging from simple dose calculations to the visualization of the whole body displaying damaged or affected organs and supplying experts with full diagnostic reports [3] . The benefits of a medical imaging examination in terms of its ability to yield an accurate diagnosis depends on the quality of both the image acquisition and the image interpretation. During the past century, radiology has grown tremendously due to advances in image detector systems and computer technology [4][5][6] . Moreover, whether computer-based learning [7] improves newly acquired knowledge and is an effective strategy for teaching prenatal ultrasound diagnostic skills to medical students when compared with instruction by traditional paper-based methods. This study presents a simulation technique to evaluate and quantify Luteinizing Hormone (LH) release from pituitary gonadotropes in response to short pulses of Gonadotropin-Releasing Hormone (GnRH). The work starts with explanatory introduction to this phenomena and then proposing an algorithm to do all the needed evaluations.
Reproductive hormones are characterized by their pulsativity. Quasi-regular pulses of Gonadotropin Releasing Hormone (GnRH) are released from the hypothalamus at variable frequency, often approximately once an hour and travel through the portal circulation to the pituitary [1,2] . They stimulate the coincident release of Luteinizing Hormone (LH). Normally, LH release is also quasi-regular and pulsatile, with amplitude that varies with the input.
However, high-frequency or continuous administration of GnRH disrupts the LH-GnRH relationship and LH production is reduced, a process termed desensitization [3] . Reduced LH retards follicular development, steroid production and ovulation. GnRH and LH also display different pulse shapes. Typically, in vivo, GnRH input pulses closely resemble square waves, the LH secretion profile exhibits an initial peak followed by an exponential decay [1] .
Given the complexity of the GnRH-LH relationship, it is useful to augment experimental study with a theoretical approach, formulating and fitting a model. The model developed here organizes the relationships among variables, describing system behavior concisely and quantitatively. Oscillations are not included in the model with the assumption of having one pool of realizable LH. This computational model explains the qualitative features of LH release in response to GnRH pulses having various durations and different concentrations both in the presence and absence of external Ca 2+ . It also provides estimates for unobserved system variables over time and can be used to predict future system response. By formulating the model in a quantitative framework, statistical procedures can be used to assess agreement between the model results and the experimental data, testing specific physiological hypotheses. Here these tests provide evidence of receptor desensitization.
It is known that Luteinizing Hormone (LH) and Follicle-Stimulating Hormone (FSH) secretion by gonadotropes located in the anterior pituitary is stimulated by Gonadotropin-Releasing Hormone (GnRH), a decapeptide that is released by the hypothalamus. Conn et al. [4,5] suggested that dimerization of the GnRH receptors on the surface of the gonadotropes was sufficient to initiate the release of LH and it has been established that this event occurs in response to agonist occupancy of the receptor [6] . Although, GnRH induces oscillations in CAC in gonadotropes Via Voltage-Gated Ca 2+ Channels (VGCCs) located in the plasma membrane [7] , its initial rise in response to GnRH is independent of extracellular Ca 2+ but is due to the rapid release of Ca 2+ from the ER. For low and medium concentrations of GnRH, the initial spike and the subsequent plateau level of CAC is increased and in the frequency of the subsequent oscillations with increasing concentrations of GnRH, however, as concentration of GnRH increases, the frequency increases with no further increase in the initial spike level of CAC. No develop model as yet applies to the CAC oscillations of gonadotropes [8] .
Although arachidonic acid is also involved in the signaling system by which GnRH causes LH release [9] , it is not included in the present model because an inhibitor of diacylglycerol lipase causes dose-dependent inhibition of LH release, without affecting the ability of arachidonic acid to facilitate LH release [10] .
Multiple proteins are also involved in GnRH signaling [11] . G protein is involved in the initiation of the GnRH-activated signaling pathway is inhibited by cholera toxin, but not by pertussis toxin and is, therefore, a G s guanyl nucleotide-binding protein [12] . Moreover, obtained data indicated two pools of LH in gonadotropes and redistribution of LH by GnRH from a non-releasable pool to the releasable pool.
Mathematical models have been developed to explore some of the Effect of dimerization on the response of the gonadotropes to various concentrations of GnRH and some related peptides were explored using some mathematical models [13,14] . However, such models did not include the interaction of the dimerized receptors with G proteins or the subsequent complex intracellular signaling systems but they mainly focused on the kinetics of receptor binding and dimerization. With increasing understandings of the signaling systems between GnRH binding and LH release [7,15] . perfusion systems development and optical methods to study the changes in cytosolic Ca 2+ content of individual gonadotropes has allowed much new data to be obtained on the changes in cytosolic Ca 2+ concentration (CAC) and the rates of release of LH in response to short pulses of GnRH.

MATERIALS AND METHODS
The Hypothesis: Response to short pulses of GnRH is only considered in order to avoid the effects of changes in gene expression that are known to occur in response to long exposure to GnRH [16] . This has made it possible for the proposed model to generate pituitary cell lines expressing different concentrations of GnRH receptors [17] , thus, to obtain data on the effect of receptor number on rates of release of LH and of the change in cytosolic concentrations of IP 3 , the signaling compound that activates release of Ca 2+ from the Endoplasmic Reticulum (ER) [18] .
The proposed model includes the interaction of the dimerized receptors with a G protein in the cell membrane, the consequent release of IP 3 and, therefore, of Ca 2+ from the Ca 2+ stores in the ER (CAER), the opening of Ca 2+ channels in the plasma membrane and the subsequent active transport of Ca 2+ back into the ER and into the extracellular fluid. Many features of the complex pathway between the binding of GnRH to its receptors and LH release were omitted. Furthermore, it focuses on data concerning LH release in response to relatively short pulses of GnRH and the effects of varying GnRH receptor number (GnRHR).
The CAC (and the CAER) were modeled as smooth functions of time that mimic the initial spike and the subsequent average Ca 2+ levels. Although the rate of CAC oscillations increases with increasing IP 3 concentration, the IP 3 level, on which the GnRHinduced Ca 2+ responses depend, does not oscillate [19] . Besides, low CAC facilitates and high CAC inhibits the release of Ca 2+ from the ER were evident. Also, it is chosen to treat LH release as occurring from a single pool because the kinetics of transfer of LH from the non releasable pool to the releasable pool has not been studied.
As the estimated for the number of GnRHRs per gonadotrope is 10 4 , then assuming that there are 10 5 -10 6 cells/ml in a typical experiment on isolated gonadotropes, the GnRHR concentration would be 1.5-15 pM. Therefore R 0 = 0.01 nM is chosen as standard total receptor concentration. Although, the G proteins per gonadotrope concentration is not known but, it is ~10-fold or more as compared to the number of receptors [20] , It is also chosen that total concentration of G protein to be GQ 0 = 0.1 nM.
Formalization: The mathematical model proposed here for the stimulation by GnRH of LH release by pituitary cells can be divided into three phases. The first phase treats the hormone bindings to the receptor, formation of dimmers, the production of effectors and IP 3 occurrences. In the second phase, IP 3 -regulated channels on the ER allow Ca 2+ to be released and subsequently pumped back into the ER. Then in the third phase, the voltage-sensitive cell membrane Ca 2+ channel, leakage Ca 2+ channels, the cell membrane Ca 2+ pump and the release of LH are described. The three stages of calculations are outlined in the following.
First Phase: Hormone binding to IP 3 formation: The concentration of GnRH in the surrounding medium as a function of time is denote by H(t) , where time t is measured in minutes. The receptor concentration per unit volume will be denoted by R. It is assumed that GnRH binds to the receptors via a simple reversible reaction with certain binding and unbinding rates k 1 and k -1 , respectively. The resulting hormone-receptor bound complex HR interacts with each other reversibly forming dimmers, whose concentration is denoted by HRRH with binding and unbinding rates k 2 and k -2 , respectively. These dimmers react also reversibly with G protein, denoted by GQ, producing effectors with concentration E at the rate of k 3 (and k-3 for reverse operation). This effecter E is assumed to represent phospholipase C. Therefore, looking at all the above reversible reactions together into an inter-related equation, one will end up with the following equation: This equation is implemented in the diagram of Fig. 1 and the concentration changing rates for R, H, HR, HRRH, GQ and E values are computed by partial differentiation with time. In order to produce inositol 1, 4, 5-triphosphate (IP 3 ). This is simply achieved by assuming that IP 3 is proportional to the effectors E concentrations. The subsequent metabolism of IP 3 in gonadotropes is complex [18] with unknown kinetics. However, it can be assumed simply that IP 3 is converted to inactive metabolites at a rate proportional to its concentration, thus yielding the following differential equation: 5 5 IP3 The total receptor concentration R 0 , or the number of GnRHRs per gonadotrope can be taken = 0.01 nM [17,22] and number of G proteins per granadotrope to be GQ = 0.1 nM [20] .  Table 1. Degradation of GnRH has made it difficult to measure the affinity constant for the binding of the naturally occurring ligand to its receptors, therefore no precise value is unavailable and hence a value of ~0.7×10 9 M −1 has been suggested [22] .
The magnitudes of k 1 and k −1 were chosen so that most of the binding occurs within 30 sec and the affinity constant is 0.5×10 9 M −1 . The values of the rate constants for dimerization k 2 and k -2 were chosen such that there is a high and rapid tendency to dimerize. Also, the rate constants k 3 and k -3 were chosen so that the binding of the dimer to the G protein is rapid and has high affinity. The rate constants k 5 and k -5 were chosen so that IP 3 approaches its steadystate level in response to H = 0.1 nM in ~2.5 min and to H = 10 nM in ~0.5 min, consistent with the short-pulse data of Morgan et al. [18] .
Applying a 5 min pulse of hormone concentrations, GnRH H = 0.1, 1 and 10 nM, the responses of the variables R, HR, HRRH, E and IP 3 are plotted in Fig. 2. At H = 0.1 nM, ~95% of the receptors R A are unoccupied and thus the amounts of of concentrations are all very small, see HRRH A , E A and IP 3A . As H is increased to 1 nM, about two-thirds of the receptors are occupied and much more is produced, see HRRH B , E B and IP 3B . With a further 10-fold increase in H, ~96% of the receptors R are occupied and further significant increases occur in other concentrations, see HRRH C , E C and IP 3C . At this high hormone concentration, the formation of HR is so rapid that HR first peaks (at ~0.1 min) and then declines as HRRH is formed. The concentrations approach their steady-state values in~1 min and then relax back to their initial values within 2-5 min after the hormone is removed at t = 5 min for all the three cases.

Second phase: Release of Ca 2+ from endoplasmic reticulum ER:
When inositol 1,4,5-trisphosphate, IP 3 binds to receptors on the Endoplasmic Reticulum (ER) membrane, Ca 2+ is stored in the ER and released. The detailed dynamics of these receptors play a major role in the cytosolic Ca 2+ oscillations and lots of research is carried out in this context [23] . However, an assumption is adopted to give a time course to fraction of open channel CHO instead of modeling these dynamics [24] : This fraction of open Ca 2+ channels depends on IP 3 concentration in a Michaelis-Menten-type saturating fashion. The fraction in the first factor on the right approaches a maximum of 1 for high IP 3 concentrations. The factor βte (1-βt) approaches a maximum of 1 when t = 1/β. The parameters α and β are chosen as 2 nM -1 and 4 min −1 , respectively [25] . The maximum probability of opening Ca 2+ channels or the factor CHO reaches its maximum of 0.6 at 0.25 min, which is in full consistent with the data of Ramos-Franco et al. [26] . Equation 3 is plotted; see It is assumed that the release of Ca 2+ is proportional to fraction of open channel CHO, difference between ER Ca 2+ concentration (CAER) and Cytosolic Ca 2+ concentrations (CAC) and the some intrinsic rate constant kk 6 [24] . The ERR equals to the rate constant k 6 plus some complicated and nonlinear function related to CAC. Also, Ca 2+ is pumped back into the ER at a rate jointly proportional to CAC and to difference between resting concentration of Ca 2+ in the ER (ERUL) and CAER. k -1 = 5 min -1 , k 2 = 2,500 nM -1 · min -1 , k -2 = 5 min -1 , k 3 = 4,000 nM -1 · min -1 , k -3 = 200 min -1 , k 5 = 2 × 10 7 min 1 , k -5 = 10 min -1 k 6 = 1 µM -1 · min 1 , k -6 = 5.0 min -1 , k 66 = 10 µM -1 · min -1 , k 666 = 0, k 7 = 2.2 µM/min, k 8 = 0.4 nM -1 · min -1, k 88 = 0, k 888 = 0, k 66 = k 88 CAC -k 888 CAC 2 k 9 = 0.0002 min 1 , k 10 = 5 ng/min Therefore, the rate of ER Ca 2+ concentration can be calculated first then the estimated dynamics rate of Cystolic Ca 2+ Concentration in the absence of other mechanisms would be calculated, as summarized in the block diagram of Fig. 4. No experimental measurements of the ER volume ratio to any cell volume, but 1/20 is found reasonable. From [8,24,27] , the resting level of cytosolic Ca 2+ is 0.05-0.2 µM [8] , the rate constants (in the absence of H) is chosen as 0.1 µM in the absence of hormones, ERUL = 40 µM, CAER = 20 µM. So by choosing k 6 = 5 µM −1 min −1 and suitable other factors, we would obtain an initial rate of decrease of CAER of about 80 µM min −1 . Then, if all channels were open or CHO = 1, the initial rate of increase of CAC would be 1/20th of that, i.e., 4 µM min −1 . However, these are overestimates, because, in fact, for the first few seconds only a small fraction of the ER channels are open and the maximum probability of opening is 0.6, this means that the initial spike in CAC occurs in<1 min and reaches a maximum value at 0.2-1 µM. these results are shown in Fig. 5 and they are in agreement with experimental data. Moreover, initial spike of CAC can be attributed to the highly weak pumping of Ca 2+ back into the ER as compared with the release of Ca 2+ from the ER.
The LH release rate is plotted for a 5 min pulse of hormone of three different concentrations, 0.1, 1 and 10 nM in Fig. 5. When H = 0.1 nM, an initial rise from 0.1 to~0.19 µM is noticed followed by a gradual decrease, but when H = 1 nM, the initial spike occurs more rapidly and goes to a much higher level, attaining a plateau while GnRH is still present. However, at H = 10 nM, the initial rate, the peak level and the plateau level are somewhat higher, but the qualitative behavior remains the same. On termination of the pulse at 5 min, CAC returns to its resting level within ~5 min. In all three cases, the LH production rate follows the CAC.   [28] that the rate is facilitated by low concentrations of CAC and inhibited by high concentrations. This dictates that Ca 2+ pumps in the cell membrane obey second-order Michaelis-Menten kinetics (rate constant k 7 ) and that Ca 2+ leakage from outside to inside (rate constant k 9 ) is a simple first-order process. The addition of these Ca 2+ fluxes the differential equation and finally using second-order Michaelis-Menten kinetics, rate of release of LH is computed depends on CAC, as shown in the diagram of Fig. 6.

RESULTS
In calculating the CAC rate, at the resting level of CAE = 1 mM (normal range of plasma free Ca 2+ concentration), a balance between k 7 and k 9 , so that the rate of pumping out of the cell equals the rate of leakage into the cell . The second-order kinetics may be chosen in order that, in the absence of GnRH, the low CAC fluctuations about resting level and cause low baseline LH release. Obviously, CAC returns to its resting state within a few minutes after the removal of GnRH. Value of k 8 is chosen such that typically observed elevated CAC values after the initial spike while GnRH is still present, while k 88 and k 888 are set to 0 and k 10 is chosen = 5 ng min −1 in order to give consistent LH release with observed rates. The release rate of LH is plotted against the time for the 5 min of GnRH ulse in Fig. 7. Experimental results using 7 min pulse [8]

DISCUSSION
The observed results listed in Fig. 7 for the release rate of luteinizing hormone LH due to a pulse for certain duration gonadotropine-releasing hormone (GnRH) are in full agreements with the studies of Stoljilkovic et al. [8] . In the presence of Ca 2+ , an initial rapid rise in LH release rate is followed by a rapid decline, a brief pause and then decay to zero in few minutes. They are also in consistent with the experimental results of Chang et al. [29] . Although the exposure times to different concentration of GnRH for these experiments were different, the main result profiles were similar. In the absence of Ca 2+ , the peak level of LH release is considerably low and there is no pause, typically as observed by Chang et al. In both cases, the rate of LH release is closely correlated to cytosolic Ca 2+ .
Different gonadotropin-releasing hormone GnRH concentration levels were implemented in Stoljilkovic et al. [8] experiments for 7 min pulse duration. For each GnRH concentration, the rate increases rapidly, reaches a peak, then declines rapidly to a plateau level during the pulse and then declines to resting level within a few minutes. The rate of LH release tracks the Cytosolic Ca 2+ concentration as shown in Fig. 7b, these results are in full consistence with results shown in Fig. 7a.
It must be noted that in the proposed study here, some important short-term mechanisms were ignored.

CONCLUSION
The proposed model demonstrates that the major features of the short-term responses can be simulated and fairly understood using only the mechanisms that are included in the model, such as receptor binding, dimerization, interaction of the dimerized receptor with a G protein, production of an effector that opens voltage sensitive channels in cell membrane and catalyzes the formation of IP 3 , which opens the Ca 2+ channels in the ER and the Ca 2+ -dependent release of LH.
This computer aided simulation and evaluation model also addresses many of the significant short-term features such as binding of GnRH to its receptors and release of LH (and FSH) from gonadotropes including the shape and time course of LH release in response to GnRH pulses which turned out to be in full agreement with seen experimentally, in the presence and absence of external Ca 2+ . Prediction of the effects of varying concentrations of Ca 2+ concentration or receptor number and comparison with experimental data was made possible. Further elaboration in feature visualization using efficient computer imaging techniques can be thought of as future extension for the current study.