Professor Ray K Iles
email: r.iles@mdx.ac.uk

P J Delves

S A Butler

Does hCG or hCG play a role in cancer cell biology?

Keywords: hCG, hCGß, differentiation, oncogenesis, apoptosis, cancer

A c c e p t e d M a n u s c r i p t 2

Summary

The role that hCG might play in the oncogenic process in cancer is certainly complex.

We know that the expression of hCG and its beta subunit is a widespread phenomenon which has been described in many cancer subtypes. However, hCG's involvement in breast cancer has been antithetical: the detection of ectopically expressed hCG() by breast tumors has been employed as a biomarker of malignancy, and hCG has been proposed as a ligand vehicle for toxic drugs, with the aim of targeting the LH/hCG receptor which is reported to be expressed by malignant breast tissue. However, it has also been proposed that hCG is a protective agent against the development of breast cancer, leading some to advocate hCG administration to non-pregnant women as a prophylactic measure against cancer. Nevertheless, suggestions that hCG is involved in the angiogenesis, metastasis and immune escape that are central to cancer progressionare phenomena which clearly apply to breast cancer. Indeed, a tumor vaccine based upon hCG has very recently been shown to protect against mammary tumors in mice.

We propose that this apparent paradox is resolved if the free beta subunit of hCG produced by tumors acts as an autocrine anti-apoptotic and angiogenic growth factor, whilst intact heterodimeric hCG, as in pregnancy, is part of developmental signaling that initiates tissue differentiation (including breast ductal tissue development), and hence reduces the population of stem-like cells which are susceptible to oncogenic factors.

Introduction

Ectopic human chorionic gonadotropin (hCG) expression by non-gestational tumors was noted as early as 1904 (DeJewitzi reported a case in which he concluded that a gonadotropin was produced by a bladder cancerreviewed by [START_REF] Iles | Human chorionic gonadotrophin expression by bladder cancers: biological and clinical potential[END_REF] and was reported by Ascheim and Zondek following the introduction of their seminal bioassay for hCG in 1927 [START_REF] Ascheim | Hypophysenworderlappen hormon und ovarialhormon im Harn von Schwangeren[END_REF][START_REF] Ascheim | Die Schwangerschaftsdiagnose aus dem Harn durch Nachweis der Hypophysevordem Lappenhormon[END_REF][START_REF] Zondek | Versuch einer biologischen (hormonalen) Diagnostik beim malignen Hodentumor[END_REF][START_REF] Zondek | Hormone des Ovariums und des Hypophysenvorderlappens[END_REF]. Up to the late 1970s, the incidence of ectopic hCG expression by common epithelial cancers could only be gleaned from a number of case reports. However, since the 1980s publications describing the ectopic expression of hCG have grown exponentially, and such findings are now generally reported in terms of the percentage of hCG positive patients with a given type of epithelial cancer (reviewed by [START_REF] Iles | Human urothelial carcinomas -A typical disease of the aged: The clinical utility of chorionic gonadotrophin in patient management and future therapy[END_REF]. A major reason for this improved appreciation of the frequent occurrence of ectopic hCGand, in particular, hCG expression -in many cancers was the introduction by [START_REF] Tanaka | A radioimmunoassay which specifically measures human chorionic gonadotropin in the presence of human luteinizing hormone[END_REF] of a radioimmunoassay which detects the free subunit of hCG, in addition to the intact hormone. Extensive immunochemical characterization has shown that, whilst intact hCG is abundantly produced by the placenta and by germ cell tumors, it is the free  subunit (independently of the common glycoprotein hormone subunit [GPHα] that is predominantly produced by common epithelial tumors (Iles & Chard, 1989;Iles et al., 1990a, Iles et al., 1990b, Stenman et al 2004). Ectopic production of free hCGβ by bladder carcinoma is well established (reviewed by [START_REF] Iles | Human urothelial carcinomas -A typical disease of the aged: The clinical utility of chorionic gonadotrophin in patient management and future therapy[END_REF], Hotakainen et al 2002b), and its expression has also been shown in a broad range of common epithelial cancers including breast, cervical, endometrial, prostate, lung, colon, oral/facial tissue and stomach (Table 1). The incidence is very variable but on A c c e p t e d M a n u s c r i p t 4 average ~30% of epithelial cancer patients will have elevated serum levels, 37% elevated urinary levels and 48% of tumors will stain hCG positive. Intact hCG heterodimer may also be produced by various common epithelial tumors in an unpredictable fashion, most often by lung cancer and hepatoblastoma [START_REF] Stenman | Human chorionic gonadotropin in cancer[END_REF]. In bladder cancer, the expression of the hCG heterodimer is associated with carcinomatosis. However, the free subunit is the most abundant form of immunoreactive hCG expressed by non-germ cell tumors [START_REF] Iles | Human chorionic gonadotrophin expression by bladder cancers: biological and clinical potential[END_REF]. An association between hCGβexpression, and anaplastic and advanced disease, has been noted in virtually all the published studies (Table 1) and many authors have commented on the aggressive nature of the hCGβ-positive tumors (Reviewed previously, [START_REF] Delves | Designing a new generation of anti-hCG vaccines for cancer therapy[END_REF]

The biological action of hCG on epithelial tumors

The association of hCGβ expression by epithelial cancers with a tendency for the tumor to resist radiotherapy and develop metastases has, in the past, been largely dismissed as an epiphenomenon of no functional significance for oncogenesis and/or tumor spread. Thus hCG (and by association, hCGβ) expression has largely been regarded solely as a marker of the presence of pluripotent stem/germ cells. In 1996, it was demonstrated that bladder cancer cell numbers increased in a dose-dependent manner following incubation with hCGβ [START_REF] Gillott | The effects of β-human chorionic gonadotrophin on the in vitro growth of bladder cancer cell lines[END_REF]. No effect could be seen following treatment with intact hCG, GPHα or hCGβ core fragment (hCGβcf, a naturally-occurring proteolytic cleavage fragment of hCGβ), but the growth stimulation was abolished by co-incubation with anti-hCGβ antiserum [START_REF] Gillott | The effects of β-human chorionic gonadotrophin on the in vitro growth of bladder cancer cell lines[END_REF]. Furthermore, addition of these antibodies to the media used to culture bladder cancer cell lines inhibited the growth of the 5637, ScaBER and RT112

A c c e p t e d M a n u s c r i p t 5 cell lines that themselves produce free hCGβ, but did not affect the growth of the T24 bladder cell line that does not produce hCGβ [START_REF] Gillott | The effects of β-human chorionic gonadotrophin on the in vitro growth of bladder cancer cell lines[END_REF]. It was also noted that the cell lines least affected by the β-subunit were those that secreted the highest concentrations of hCGβ, suggestive of autocrine stimulation. The hCGβ/hLHβ gene cluster has been shown not to be amplified or rearranged in bladder tumor cells, indicating that ectopic hCG expression is more likely to be the result of altered gene regulation (Iles et al., 1989). The inhibitory effects of anti-hCG antibodies on the growth of tumor cells has also been demonstrated in vivo in athymic (nude) mice with choriocarcinoma [START_REF] Bagshawe | A cytotoxic agent can be generated selectively at cancer sites[END_REF]. Furthermore, anti-sense oligonucleotides directed to the hCG β-chain have antiproliferative effects on prostate and cervical tumor cell lines in vitro [START_REF] Devi | Inhibition of human chorionic gonadotropin beta-subunit modulates the mitogenic effect of c-myc in human prostate cancer cells[END_REF][START_REF] Jankowska | Reduction of human chorionic gonadotropin beta subunit expression by modified U1 snRNA caused apoptosis in cervical cancer cells[END_REF].

Other studies showed that the increase in the cancer cell population in response to hCG was brought about by inhibition of apoptosis rather than by stimulation of cell replication [START_REF] Butler | The increase in bladder carcinoma cell population induced by the free  subunit of human chorionic gonadotrophin is a result of an anti-apoptosis effect and not cell proliferation[END_REF]. Furthermore, given that the LH/hCG receptor is not expressed by these tissues, and that only intact hCG can stimulate the LH/hCG receptor, any activity observed must be mediated by an as yet unidentified pathway [START_REF] Davies | Ectopic expression of glycoprotein hormones and their receptors by urogenital cancers[END_REF].

hCG is a member of the cystine knot growth factor/TGF superfamily

In 1994 the X-ray crystallographic structure of hCG was solved [START_REF] Lapthorn | Crystal structure of human chorionic gonadotrophin[END_REF][START_REF] Wu | Structure of human chorionic gonadotropin at 2.6 Å resolution from MAD analysis of the selenomethionyl protein[END_REF]. The most striking feature was the arrangement of three disulfide bridges in the centre of each subunit. The positions of the three cystines are almost identical in both subunits, where two disulfides bridge the anti-parallel strands of the peptide chain A c c e p t e d M a n u s c r i p t 6 forming a central loop, through which the third disulfide passes. This structure had been previously identified in a number of factors which are consequently designated the 'cystine knot growth factor' (CKGF) family. They include transforming growth factor  (TGFβ), platelet-derived growth factor B (PDGFB) and nerve growth factor (NGF). Since the identification of the structure of hCG in 1994, the family of growth factors has expanded, and now includes the vascular endothelial growth factors (VEGFs) and placental growth factor (PlGF) -despite sharing structural similarities the various members of the CKGF family carry out quite distinct functions (see figure 1). These additions to the CKGF family have further increased the intriguing possibilities for crosstalk between multiple growth regulatory systems, and ectopic expression of molecules containing a cystine knot may give rise to disruptions that favor cancer development. In fact, homology modeling using current data on the structure of hCG as obtained from crystal structure [START_REF] Lapthorn | Crystal structure of human chorionic gonadotrophin[END_REF] reveals quite striking similarities to some CKGFs in particular VEGF and TGFβ (see figure 2).

Cross talk of hCGβ with CKGF receptors

The evolution of the hCG/LH gene cluster is fascinating, and the discovery that it is interspersed by another series of pseudogenes of the cystine knot growth factor familyneurotropin -further points to its origin from, and potential to interact with, other cystine knot growth and endocrine systems [START_REF] Rull | Expression of beta-subunit of HCG genes during normal and failed pregnancy[END_REF] (see Figure 3) Whilst the hCG heterodimer binds with high affinity (K d 0.1-1.0 nM, [START_REF] Xie | Extracellular domain of lutropin/choriogonadotropin receptor expressed in transfected cells binds choriogonadotropin with high affinity[END_REF] to the CG/LH receptor, the free  subunit essentially does not stimulate the receptor. It has A c c e p t e d M a n u s c r i p t 7 been demonstrated that the free subunits will also bind the LH receptor, albeit with a much lower affinity [START_REF] Pierce | Glycoprotein hormones: structure and function[END_REF] and some studies have shown that upon receptor binding, the α subunit of hCG is internalized, and is associated with stimulation of adenylate cyclase activity. This implies that the binding model is not simple, and that hormone/receptor binding involves a topographical determinant extending across both ligand subunits [START_REF] Ryan | The glycoprotein hormones: recent studies of structure-function relationships[END_REF]. If intact hCG is required for LHR binding and stimulation, and LHR is not found in cell lines shown to respond to free hCGβ, then the biological activity described above for the free  subunit cannot be accounted for by signaling through this receptor. However, the structural homology with the cystine knot growth factors suggests that free hCG cross-interaction with other cystine knot growth factor receptors may occur.

One early feature of malignancy is neovascularisation, allowing the tumor to grow to a mass greater than 2mm. The process of vascularisation is regulated in part by cystine knot growth factors, in particular vascular endothelial growth factors (VEGFs). The development of blood vessels by the early and late placental tissue is distinct and involves VEGF receptor expression and stimulation [START_REF] Demir | Sequential expression of VEGF and its receptors in human placental villi during very early pregnancy: differences between placental vasculogenesis and angiogenesis[END_REF]. Recent studies have suggested an interplay between placental VEGF and hCG in formation of placental vasularture [START_REF] Herr | HCG in the regulation of placental angiogenesis. Results of an in vitro study[END_REF], Reisinger et al., 2007). Indeed, hCG itself has been suggested to be a novel angiogenic factor [START_REF] Zygmunt | Characterization of human chorionic gonadotropin as a novel angiogenic factor[END_REF] and we linked immunohistochemically detected hCG expression by cervical cancer with the extent of tumor vascularisation [START_REF] Li | hCG beta expression by cervical squamous carcinoma--in vivo histological association with tumour invasion and apoptosis[END_REF]. In addition, [START_REF] Arrieta | Serum human chorionic gonadotropin is associated with angiogenesis in germ cell testicular tumors[END_REF] have recently linked serum hCG levels to neo-vascularisation of non-seminomatous testicular germ cell tumors. Therefore, Given the structural homology between hCG  and the cystine knot TGFβ is encoded by a locus adjacent to that of the hCGβ-hLHβ gene cluster on chromosome 19 [START_REF] Fujii | Transforming growth factor β gene maps to human chromosome 19 long arm and to mouse chromosome 7[END_REF], and high levels of TGFβ and its receptors are expressed by urothelial tumors, despite the growth inhibitory effects of this cytokine [START_REF] Coombs | Reduced expression of TGF is associated with advanced disease in transitional cell carcinoma[END_REF][START_REF] Eder | Transforming growth factors-1 and 2 in serum and urine from patients with bladder carcinoma[END_REF][START_REF] Eder | Expression of transforming growth factors-1, 2 and 3 in human bladder carcinomas[END_REF]. This has long been a puzzle to cancer researchers; disruption of the TGFβ apoptosis signaling pathway has been considered an important feature of oncogenesis, i.e. escape from a fundamental mechanism of regulating the normal tissue epithelial cell population [START_REF] Barcellos-Hoff | Transforming growth factor  and breast cancer: Mammary gland development[END_REF]. We have previously shown that exogenous addition of hCGβ to bladder cancer cell cultures reverses the apoptosis-inducing effects of TGFβ [START_REF] Butler | The increase in bladder carcinoma cell population induced by the free  subunit of human chorionic gonadotrophin is a result of an anti-apoptosis effect and not cell proliferation[END_REF]. Ectopic expression of free hCGβ may thus be acting as an autocrine/paracrine growth-stimulating factor by virtue of its ability to block TGFβ's action although we have never published receptor binding data. Its structural homology to TGFβ suggests that free hCGβ may bind to a component of the TGFβ receptor complex (for example, TGFRII) or affect other downstream apoptopic signals.

Therapeutic targeting and selection of hCGβ expressing tumors

The molecular pathway proposed above is a reversal of the Her2 -Herceptin® molecular mechanism, but the parallels are somewhat compelling. The step that propelled Her2 from being a marker of poor prognosis was that it acted as both the selection criteria for a A c c e p t e d M a n u s c r i p t 9 therapy, and was itself the target of a novel humanized antibody treatment. To date, no infusion of antibody has been attempted to treat hCGβ-expressing tumors, but an anti-hCGβ carboxy terminal peptide (CTP) vaccine -a much greater step -has completed a phase II clinical trial in patients with colorectal cancer [START_REF] Moulton | Active specific immunotherapy with a beta-human chorionic gonadotropin peptide vaccine in patients with metastatic colorectal cancer: antibody response is associated with improved survival[END_REF]. A high anti-hCG titre response to the vaccine was associated with improved survival. However, several studies have indicated that only 15-41% of patients with colorectal cancer have increased levels of circulating hCGβ [START_REF] Stenman | Human chorionic gonadotropin in cancer[END_REF][START_REF] Moulton | Active specific immunotherapy with a beta-human chorionic gonadotropin peptide vaccine in patients with metastatic colorectal cancer: antibody response is associated with improved survival[END_REF], and the vaccine was administered to patients irrespective of whether they were 'hCG positive' or not. A vaccine which stimulates immune responses against hCGβ would be expected to benefit only those patients in whom hCGβ is aberrantly expressed, suggesting that the vaccine might have benefited a higher proportion of the trial group if a lack of circulating hCG had been applied as an exclusion criterion for entry into the trial. Indeed, the Her2 paradigm clearly showed that Herceptin would never have displayed the small but significant advantage in breast cancer patient survival if it had been given to all breast cancer patients irrespective of their Her2 expression status [START_REF] Baselga | Clinical trials of Herceptin(trastuzumab)[END_REF].

Taking many of these lessons on board Celldex Therapeutics (Needham, MA, USA) have developed an anti-hCGβ cancer vaccine based on their patented technology for targeting antigen presenting cells (APC) - [START_REF] Keler | Antibody-targeted vaccines[END_REF]. Termed CDX-1307,

An hCG-specific anti-tumor vaccine engineered to optimize the immune response

The hCG vaccine administered in the phase II clinical trial in patients with colorectal cancer [START_REF] Moulton | Active specific immunotherapy with a beta-human chorionic gonadotropin peptide vaccine in patients with metastatic colorectal cancer: antibody response is associated with improved survival[END_REF] consisted of the CTP of hCG conjugated to diphtheria toxoid. Although -unlike the rest of the hCG -chain -the CTP has the advantage of lacking sequence homology with the other members of the glycoprotein hormone family, it is known that it gives rise to relatively low affinity antibodies [START_REF] Talwar | Vaccines and passive immunological approaches for the control of fertility and hormone-dependent cancers[END_REF].

Furthermore, it clearly will possess a restricted range of hCG-specific epitopes compared to the full length -chain, which contains a total of five unique epitopes: the linear epitopes β8 and β9 located within the CTP and the discontinuous epitopes β1, β6 and β7 found in the rest of the hCG β-chain [START_REF] Berger | The ISOBM TD-7 Workshop on hCG and related molecules. Towards user-oriented standardization of pregnancy and tumor diagnosis: assignment of epitopes to the three-dimensional structure of diagnostically and commercially relevant monoclonal antibodies directed against human chorionic gonadotropin and derivatives[END_REF]. These epitopes are not present on other members (LH, FSH and TSH) of the glycoprotein hormone family. The other epitopes on the hCG -chain are designated β2, β3, β4 and β5 and are cross-reactive with similar epitopes on the LH -chain. Using oligonucleotide-directed site-specific mutagenesis, we made a number of amino acid replacements in the hCG β-chain with the aim of eliminating the cross-reactive epitopes [START_REF] Jackson | Identification and selective destruction of shared epitopes in human chorionic gonadotropin beta subunit[END_REF][START_REF] Delves | Designing a new generation of anti-hCG vaccines for cancer therapy[END_REF].

One of the mutants, with only a single amino acid replacement (a glutamic acid substituted for the arginine at position 68), dramatically refocuses the immune response away from the immunodominant cross-reactive epitope region towards to hCG-specific 1, 6 7 and CTP epitopes, with very little cross-reactive antibody production in both A c c e p t e d M a n u s c r i p t 11 rabbit and mouse immunogenicity studies [START_REF] Porakishvili | Elimination of luteinizing hormone cross-reactive epitopes from human chorionic gonadotropin[END_REF][START_REF] Chiesa | Refocusing of B-cell responses following a single amino acid substitution in an antigen[END_REF].

Strikingly, the antibodies produced in response to this mutant (hCGR68E) bind avidly to the CTP [START_REF] Porakishvili | Elimination of luteinizing hormone cross-reactive epitopes from human chorionic gonadotropin[END_REF]. This molecule is therefore a candidate for an improved hCG-specific cancer vaccine.

The paradox

Ectopic expression of hCGβ serves as a molecular biomarker identifying a subpopulation of patient tumors with poor prognosis for which molecular mechanisms have been identified in vitro and which provides the basis of a therapeutic target for antibody therapy; similar to the Her2 Paradigm. Although this paradigm may hold for epithelial tumors such as bladder, cervix and lung, it is not as clear-cut for breast cancers. However, within the current scientific literature controversy exists over the role of the pregnancy hormone hCG in the pathogenesis of breast cancer. Whilst some have demonstrated that ectopic expression of hCGβ by the breast tumor itself is a marker of poor prognosis, others have proposed hCG to be a protective agent against the subsequent development of breast cancer.

Ectopic hCGβ expression by breast cancer as a marker of poor prognosis

Studies of ectopic hCG expression by breast cancer have been intriguing (Tables 2 &3).

Most have concentrated on immunohistochemical and mRNA expression as the method of detecting ectopic expression of hCG, specifically the beta-subunit. The incidence, combining all studies of immunohistochemical detection of ectopic expression by breast A c c e p t e d M a n u s c r i p t 12 tumors, is low -24% (148/621) -compared to the general incidence of ectopic expression by epithelial cancer described in table 1. Two of the five studies have specifically stated IHC detection of hCG expression by breast tumors does not to correlate with patient prognosis, while one eluded to prognosis as expression correlated with a specific histological type of breast cancer (apocrine) [START_REF] Lee | The immunoexpressions and prognostic significance of inhibin alpha and Beta human chorionic gonadotrophins (HCG) in breast carcinomas[END_REF].

Messenger RNA studies concentrated on looking for shed blood-circulating tumor cells, and used detection of hCG mRNA as one of the markers. Combining these studies, 36-38% of patients had circulating breast cancer cells shed from the primary tumor detected as expressing hCG mRNA (Table3). However, Hotakanien et al (2002b and2007) have also suggested that false positives may arise due to poor primer design and result from LHβ expression by normal blood cells.

hCG as an anti-breast cancer therapy

The case for hCG being an agent with which to prevent cancer is encapsulated in José and Irma Russos' 2004 book, the Molecular Basis of Breast Cancer: Prevention and Treatment. The incidence of breast cancer is still rising, is the most prevalent cancer in Europe and USA, and epidemiology is increasingly indicating that many of the key risk factors are related to life-styles associated with Western cultures; refined food diet, early menarche, obesity, and that women have their first full-term pregnancies when they are older. Breast cancer is not common in some Far Eastern countries, although women from those societies are at greater risk if they migrate to the West than remain in their native country. Identifying how these risk factors impact on breast development at a molecular A c c e p t e d M a n u s c r i p t 13 level could enable preventative strategies to be established. Russo and Russo focused on the observation that early pregnancy protects against breast cancer. This is not surprising as the breast epithelia dramatically alters during pregnancy, from a stem cell rich and dormant layer to a highly differentiated multiple-structured tissue capable of synthesizing and releasing large quantities of milk in response to key specific stimuli.

They go on to propose that a treatment to artificially induce differentiation within mammary epithelial cells could represent a radical new approach to prevent much of the disease. The epidemiological data from which this hypothesis arises is supported strongly by experimental observations that a state of pregnancy -whether induced naturally, by estrogen/progesterone treatment, or by human chorionic gonadotropin -has a dramatic inhibitory effect on carcinogen-induced mammary cancer in rats. The authors argue a case for ectopic administration of hCG as a possible preventative strategy. However, although this might be a viable strategy for eliminating tumors in rats, it might be premature to consider it as preventative in humans. The contentious issue in extrapolating the result from experimental studies on rats to humans is that the rat does not possess an equivalent hCG gene and their conceptae do not produce any gonadotropic-like hormone.

It is likely that the cell-differentiating properties of sustained exposure to ovarian steroids, as occurs in pregnancy, are part of a complex and coordinated temporal staging of pregnancy-associated cytokine signaling that induces breast epithelial stem cell differentiation towards more terminally developed lactational epithelia. The further along this pathway that differentiation proceeds, the greater the reduction in potential to develop breast cancer. More recent studies suggest that breast tissue expresses LH/hCG receptors and that hCG has the direct affect of human breast epithelial differentiation (see A c c e p t e d M a n u s c r i p t 14 table 2). Indeed, for pregnancy-associated holo-hCG to have this effect, the breast epithelial tissue must express the LH/hCG receptor at functionally responsive levels.

Despite some reservations, phase I and II trials have been conducted on human patients with breast cancer receiving recombinant hCG (rhCG), and in December 2008 a phase III trial of patients at high risk for breast cancer receiving 30 doses of hCG was underway in New York (Prevention of Breast Cancer, randomized, single blind clinical trial (USA) -NCT00808522). Initial results should be published in 2010.

hCG and breast cancer

Interestingly several studies have used hCG as targeting ligand to which a toxin is conjugated [START_REF] Bodek | A novel approach of targeted ablation of mammary carcinoma cells through luteinizing hormone receptors using Hecate-CGbeta conjugate[END_REF][START_REF] Hansel | Destruction of breast cancers and their metastases by lytic peptide conjugates in vitro and in vivo[END_REF]. This suggests that such tumors express/over express LH/hCG receptors to which they respond. These two opposing positions raise the issue of whether exposure to exogenous hCG in the stimulated cycle protocol leading up to IVF have increased or decreased the risk of breast cancer.

Reviewing the literature in a meta analysis of 60,050 women undergoing IVF cycles, Salhab and colleagues found that there was no change in relative risk (rr=1.03) for the development of breast cancer [START_REF] Salhab | In vitro fertilization and breast cancer risk: a review[END_REF] . Transgenic mice over-expressing an engineered LH which contains an hCG-like carboxy terminal peptide do have an increased incidence of breast cancer, but this is probably secondary to gross hyperstimulation of ovarian function leading to anovulation, infertility, increased steroid production and ovarian granulosa cell tumorigenesis [START_REF] Mann | Consequences of elevated luteinizing hormone on diverse physiological systems: use of the LHbetaCTP transgenic mouse as a model of ovarian hyperstimulation-induced pathophysiology[END_REF]. In a similar vein of investigation, [START_REF] Chen | Do Placental Genes Affect Maternal Breast Cancer? Association between Offspring's CGB5 and CSH1 Gene Variants and Maternal Breast Cancer Risk[END_REF] examined polymorphisms around CGB gene 5 of the CG/LH gene cluster, and human placental lactogen gene of the hPL-GH A c c e p t e d M a n u s c r i p t 15 gene cluster, in an attempt to link either hormone to the molecular causation of the 4.5 fold protective effect of term pregnancy on the risk of developing breast cancer. As a large case control study, subjects and controls were selected from the New York Breast Cancer Family Register. None of the polymorphic markers for hPL were associated with an increased risk of developing breast cancer. However, the variant C allele of CGB5 rs726002 SNP had an elevated breast cancer risk (OR = 2.09) only in women bearing children at a young age, <26 years. Among women with a later age at childbirth, breast cancer risk did not differ.

Conclusion

Given the differences in structure and receptor target of intact hCG and free hCGβ, it is entirely feasible that the free beta subunit of hCG is produced by tumors and acts as an autocrine (anti-apoptotic) growth factor, whilst in pregnancy, intact heterodimeric hCG is part of a developmental signal that, among various systems, initiates breast ductal tissue differentiation and hence reduces the population of stem-like cells which are susceptible to oncogenic factors. The challenge is to truly discriminate these two actions as representing two distinct as yet to be defined molecular pathways. Case reports have been omitted as they do not add information regarding the frequency of ectopic hCGβ expression. Thirty-seven of 40 reports investigated the frequency of elevated levels of hCG (or hCGcf) in serum or urine and/or detection in tissue sections by immunohistochemistry (IHC). Positive detection ranges from 0% in renal cell cancer to 93% in small cell (SC) lung cancer; ~30% of epithelial cancer patients will have elevated serum levels, 37% elevated urinary levels and 48% of tumors will stain hCG positive. Twenty three of the studies specifically investigated cancer prognosis and 20 (87%) report a strong association between hCGβ detection and poor patient prognosismetastatic spread. Other studies have investigated mRNA levels such as prostate cancer [START_REF] Daja | Betahuman chorionic gonadotropin in semen: a marker for early detection of prostate cancer?[END_REF][START_REF] Span | Analysis of expression of chorionic gonadotrophin transcripts in prostate cancer by quantitative Taqman and a modified molecular beacon RT-PCR[END_REF] and renal carcinoma [START_REF] Jiang | Expression of beta-human chorionic gonadotropin genes in renal cell cancer and benign renal disease tissues[END_REF] but many such studies were compromised by high-level hCGβ mRNA detection in normal tissues. mRNA detection of hCG detected in breast cancer cells is described in Table 3.

  A c c e p t e d M a n u s c r i p t 8 growth factors(lapthorn et al., 1994), and its ability to dimerise[START_REF] Butler | The  subunit of human chorionic gonadotrophin exists as a homodimer[END_REF] cross talk with VEGF R is plausible.

  this vaccine completed Phase 1 development and is marketed as a treatment for colorectal, pancreatic, bladder, ovarian and breast cancers. CDX-1307 couples the entire beta chain of hCG to an APC specific adjuvant antibody and induces both humoral and cell mediated response to the complete molecule and not just the CTP. Phase 2 study in patients with invasive A c c e p t e d M a n u s c r i p t 10 bladder cancer began in 2010 but inclusion is dependant on demonstrating hCGβ expression by the tumors.

  Rong Y, Yankai Z, Wentao L, Hongxia Z, Jie W, Rongyue C, Taiming L, Jingjing L. 2006; Improved efficacy of DNA vaccination against breast cancer by boosting with the repeat beta-hCG C-terminal peptide carried by mycobacterial heat-shock protein HSP65. Vaccine. 24: 2575-2584. Yves A Muller, Hans W Christinger, Bruce A Keyt, Abraham M de Vos. (1997).Vascular endothelial growth factor: Crystal structure and functional mapping of the kinase domain receptor binding site. Proc. Natl. Acad. Sci. USA 9: 7192-7197

Figure 1 :

 1 Figure 1: Cystine Knot superfamily -Diagrammatic representation of the cystine knot hormones/growth factors grouped by functional family. Functional families in which hCG/free hCG may cross-react are expanded in detail.

Figure 2 :

 2 Figure 2: Comparative shapes of cystine knot dimers. Molecular models of the 3D structure of dimeric TGF, VEGF, NGF and a theoretical hCG homodimer.

Figure 3

 3 Figure 3 hCGβ-LHβ gene cluster ch19q13.32,-Diagrammatic representation of the relative position of hCGgenes 1,2,3,5,7,8 and LH interspersed with neurotropin pseudogenes-NTF4, 6, 6a and 6G.

Fabisiewicz

  et al., 2004 RT-PCR detects hCGß mRNA in peripheral blood circulating tumour cells of 21/55 breast cancer patients (stage 1-11 N0-N1) Chang et al., 2005 44/334 tumors stained positive for hCGβ more common in apocrine breast cancers (4/5 positive)

Figure

  Figure

Table and figure legends Table 1. Frequency and significance of hCGβ expression by common carcinomas.

Table 2 -Studies suggesting that hCG has a direct protective effect on breast tissue and anti-cancer activity. Three

 2 studies have examined LH/hCG receptor expression by breast tissue; one used such expression to target drug therapy and six studies have demonstrated hCG stimulation of secondary messengers and change in breast cancer cell gene expression as a result.

Table 3 -Studies reporting breast tumor ectopic expression of hCGβ.

 3 These studies have demonstrate expression as detected by IHC but the incidence (24%) is lower than the average (~48%, see table 1) and is not reported as associated with prognosis. hCG has been targeted as a marker for the detection of shed circulating breast cancer cells and is detected in 36-38% of cases.

Table 2 -

 2 Reference FindingsReview of potential of recombinant hCG as a protective drug Meduri et al., 2004 160 patient tumors were examined for LH/hCG receptor, 72% stained (strongly) positive. Positive staining was associated with increased metastases-free survival Jansseens et al., 2007 Pilot study of rhCG administered to patients with primary breast cancer. Excised treated tumors had reduced proliferation index and reduced estrogen &progesterone receptor expression. Metastatic breast cancer patients also received rhCG. Kocdor et al., 2009 MCF-10F cells treated with rhCG differentiated and prevented 17 -estradiol transformation of cells in collagen matrix cultures into solid masses.

	Hu et al., 1999	HCG/LH receptor expression much higher in MCF-7 cells than
		other breast cancer and tissue cell lines
	Pond-Tor et al.,	Treatment with hCG increases MCF-7 breast cancer cells
	2002	sensitivity to radiotherapy
	Gebaur et al., 2003 Conjugation of deoxrubicin to hCG provides as specific targeting
		modality for cytotoxic drug delivery to LH/hCG receptor
		positive breast cancer cells
	Guo et al., 2004	hCG significantly alters the gene expression array of breast
		cancer cell line MCF-7
	Carlsen et al., 2004 -IHC and ISH detection of LH/hCG receptor in all male breast
		tissue samples -intense IHC staining in 0/4 normal 5/13
		gynecomastia and 4/5 cancer samples
	Rao et al., 2004	hCG down regulates proliferation and matrigel invasion by
		MCF-7 by lowering responsiveness of NFkB and AP-1 to stimuli
	Janssens et al.,
	2004

Table 3 -

 3 ReferenceFindingsWachner et al., 1984 15/129 tumors stained positive (but only very few per field). Taback et al., 2001 HCGβ mRNA was detected in 19 of 25 primary breast cancer tissues and 22 of 65 peripheral blood samples from breast cancer patients at diagnosis. Primary breast cancer tumour tissue expressed functional hCGβ gene transcripts but not the mRNA for the alpha subunit Thus, although primary breast cancer also expressed mRNA for the LH/hCG receptor there could not be auto stimulation of this receptor as the alpha-beta heterodimer hCG could not be formed. Hu & Chow 2001 hCGβ mRNA was detected in 9/72 breast cancer patients peripheral blood cancer cells and detection correlated with stage of disease. Tanaka et al., 2000 Although hCG did not directly stimulate MCF-7 breast cancer cell growth, it stimulates expression of aromatase enzymes such that DHEA is converted to estradiol, and hence stimulates estrogen dependant MCF-7 cell growth. Span et al 2003 RT-PCR Analysis of hCGβ expression of 129 patients breast cancer tissue found high hCGβ mRNA levels in 65. High-level expression correlated with lower progression free survival, and was an independent marker to nodal status. Hu et al., 2003 Circulating cancer cells were isolated from blood samples taken from 49 patients. hCGβ mRNA expression transcripts were detected in 6 pre-operatively and 13 immediately post-operatively. Expression did not correlate with recurrence and metastases.

	Agnanatis et al.,	15/27 tumors stained strongly positive for hCGβ
	1992.

Acknowledgementwe grateful acknowledge the 3D modelling of hCGβ dimers by Saifur-Rehman and his bioinformatics comparisons with the cystine knot superfamily.

Professor Ray Iles is supported by the Lee Smith Foundation and Peter Luck-Hille Foundation References


Zoran Nikoloski
email: nikoloski@mpimp-golm.mpg.de

Patrick May

Joachim Selbig

Algebraic connectivity may explain the evolution of gene regulatory networks

Keywords: algebraic connectivity, evolution, gene regulatory networks, dynamics

Gene expression is a result of the interplay between the structure, type, kinetics, and specificity of gene regulatory interactions, whose diversity gives rise to the variety of life forms.

As the dynamic behavior of gene regulatory networks depends on their structure, here we attempt to determine structural reasons which, despite the similarities in global network properties, may explain the large differences in organismal complexity. We demonstrate that the algebraic connectivity, the smallest non-trivial eigenvalue of the Laplacian, of the directed gene regulatory networks decreases with the increase of organismal complexity, and may therefore explain the difference between the variety of analyzed regulatory networks. In addition, our results point out that, for the species considered in this study, evolution favours decreasing concentration of strategically positioned feed forward loops, so that the network as a whole can increase the specificity towards changing environments.

Moreover, contrary to the existing results, we show that the average degree, the length of the longest cascade, and the average cascade length of gene regulatory networks cannot recover the evolutionary relationships between organisms. Whereas the dynamical properties of special subnetworks are relatively well understood, there is still limited knowledge about the evolutionary reasons for the already identified design principles pertaining to these special subnetworks, underlying the global quantitative features of gene regulatory networks of different organisms. The behavior of the algebraic connectivity, which we show valid on gene regulatory networks extracted from curated databases, can serve as an additional evolutionary principle of organism-specific regulatory networks.

Introduction

Recent evidence from fully-sequenced genomes suggests that organismal complexity arises much more from the elaborate regulation of gene expression than by the genome size itself [START_REF] Levine | Transcription regulation and animal diversity[END_REF]. The ever-increasing throughput in experimental manipulation of gene activity coupled with the methods for quantitative assessment of perturbed transcriptome, proteome, and metabolome have begun to identify the effects of transcription factors, binding ligands, and post-translational modifications on regulated genes [START_REF] Papp | Genome-wide analysis of context-dependence of regulatory networks[END_REF]. In addition to structural information regarding the regulatory interactions, comprehensive understanding of the behavior of these interactions also requires specification of: (1) the type of regulation (i.e., activation or inhibition) [START_REF] Albert | The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in drosophila melanogaster[END_REF], (2) kinetics of interactions [START_REF] Ronen | Assigning numbers to the arrows: parameterizing a gene regulation network by using accurate expression kinetics[END_REF], and (3) the specificity of the interactions with respect to investigated tissue and/or stress condition [START_REF] Luscombe | Genomic analysis of regulatory network dynamics reveals large topological changes[END_REF]. The elucidation of complete network of regulatory interactions parametrized with kinetic information leading to a particular type of gene expression is, at present, still a challenging task even for the well-studied model organisms whose networks have been partially assembled for few selected processes, conditions or on the level of the entire genome [START_REF] Davidson | A genomic regulatory network for development[END_REF][START_REF] Lee | Transcriptional regulatory networks in saccharomyces cerevisiae[END_REF][START_REF] Shen-Orr | Network motifs in the transcriptional regulation network of escherichia coli[END_REF][START_REF] Zhang | Dynamic changes in subgraph preference profiles of crucial transcription factors[END_REF].

Nevertheless, recent theoretical investigations have established that the qualitative behavior of dynamic processes on complex networks is closely related to the structural network properties [START_REF] Conradi | Subnetwork analysis reveals dynamic features of complex (bio)chemical networks[END_REF][START_REF] Craciun | Understanding bistability in complex enzyme-driven reaction networks[END_REF][START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF][START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactors-i. the deficiency zero and deficiency one theorems[END_REF][START_REF] Madan Babu | Evolutionary dynamics of prokaryotic transcriptional regulatory networks[END_REF]. Consequently, the existing studies of gene regulatory networks have attempted to determine unifying design principles in order to understand and make biologically relevant conclusions solely from the network structure.

The reported results mainly fall in the following two categories: (1) identifying local structures, such as the basic building blocks of networks, which, if present at statistically significant concentrations, are termed motifs [START_REF] Milo | Network motifs: Simple building blocks of complex networks[END_REF][START_REF] Shen-Orr | Network motifs in the transcriptional regulation network of escherichia coli[END_REF], or particular subnetworks, called cascades [START_REF] Rosenfeld | Response delays and the structure of transcription networks[END_REF], and (2) developing models which can explain global structural properties related to the salient network properties (e.g., degree distribution, average path length, clustering coefficient) [START_REF] Albert | Statistical mechanics of complex networks[END_REF].

Network motifs are defined as patterns that occur more often in the real network than in randomized networks [START_REF] Milo | Network motifs: Simple building blocks of complex networks[END_REF]. Three types of recurring network motifs were found to describe most of the Escherichia coli and Saccharomyces cerevisiae transcriptional regulatory networks [START_REF] Milo | Network motifs: Simple building blocks of complex networks[END_REF][START_REF] Shen-Orr | Network motifs in the transcriptional regulation network of escherichia coli[END_REF]. Each of these network motifs was shown to exhibit quantitative differences in the levels and timing of gene expression [START_REF] Alon | Network motifs: theory and experimental approaches[END_REF]. Here, we are interested only in the concentration of a particular subnetwork, rather than its statistical significance. Another important feature of regulatory network architectures is the distribution of transcription cascade lengths [START_REF] Rosenfeld | Response delays and the structure of transcription networks[END_REF][START_REF] Shen-Orr | Network motifs in the transcriptional regulation network of escherichia coli[END_REF]. Transcription cascades are defined by a set of transcription factors that regulate each other sequentially. In a study of the architecture of the transcription network of Escherichia coli and Saccharomyces cerevisiae, it was observed that these networks have a strikingly shallow architecture, with most genes regulated by overlapping cascades of length one or two, which is different for developmental regulatory interactions of Drosophila melanogaster [START_REF] Rosenfeld | Response delays and the structure of transcription networks[END_REF]. Thus, it was consequently suggested that the response delay, governed by the length of cascades, may act as a design principle of gene regulatory networks.

With respect to the salient network properties, gene regulatory networks, like the majority of social and technological networks, have been shown to belong to the class of scale-free networks (also called complex networks) [START_REF] Albert | Scale-free networks in cell biology[END_REF]. In addition to the scale-free degree distribution, in which the number of highly connected nodes (called hubs) is small, these networks are marked with a large clustering coefficient and small average path length. Recently, it was demonstrated that precisely such combination of network properties allows ultra fast emergence of consensus on the network (where consensus of opinions is achieved only through communication of a network node with its direct neighbours) and synchronization of coupled oscillators robust to edge removal [START_REF] Olfati-Saber | Ultrafast consensus in small-world networks[END_REF][START_REF] Wang | Synchronization in scale-free dynamical networks: robustness and fragility[END_REF].

Whereas the dynamical properties of special subnetworks (e.g., cascades and network motifs) are relatively well understood, there is still limited knowledge about the evolutionary reasons for the already identified design principles pertaining to these special subnetworks, underlying the global quantitative features of gene regulatory networks of different organisms. Moreover, despite tremendous progress in the development of models (via simple principles/rules) that mimic the global properties of complex regulatory networks [START_REF] Albert | Statistical mechanics of complex networks[END_REF], the identified rules cannot explain the differences between regulatory architectures from different organisms.

Here, we use the adjacency matrix representation of directed regulatory networks which allow network comparison via the spectrum, i.e., the set of eigenvalues, of their corresponding matrices. Graph spectra have already been proposed as a systematic tool in computational biology and other fields, and have been applied in the analysis of undirected networks [START_REF] Banerjee | Graph spectra as a systematic tool in computational biology[END_REF]. Moreover, it has been observed that a spectrum of a graph in fact encodes information about many structural properties (e.g., average path length, diameter, treewidth) as well as local substructures which may lead to positing hypotheses regarding the evolution of the networks [START_REF] Chung | Complex Graphs and Networks[END_REF]. We demonstrate that the algebraic connectivity (smallest non-trivial eigenvalue) of the directed graph representation of gene regulatory networks is a network characteristic which can explain the difference between the regulatory networks of different organisms. Our findings about the algebraic connectivity suggests that information flow in gene networks of simple organisms is generally less spatially restricted that in the networks of more complex organisms. The behavior of the algebraic connectivity, which we show valid on gene regulatory networks from different organisms, is an additional hypothesis for the evolutionary principles of organism-specific regulatory networks giving rise to the variety of life forms.

Methods

A gene regulatory network comprises the transcription factors and two types of genes: (1) genes coding for transcription factors, called regulatory genes, and (2) genes regulated by transcription factors, known as target genes. In the case of auto-regulation and cascades, a regulatory gene also plays the role of a target gene. We model the structure of regulatory interactions as a directed bipartite graph G, in which the two partitions, V tf (G) and V g (G), comprise the transcription factors and the genes, respectively. Such a graph will be called TF-gene regulatory network. From a TF-gene regulatory network G, one can obtain a gene-gene regulatory network G . The graph G contains only those nodes from G which are in V g (G). The transformation from G to G is accomplished according to the following rule: A directed edge is drawn from node u to node v in G if and only if there is a directed path of length two between the corresponding nodes in G. The effects of this transformation on the properties of the resulting network are illustrated and disccused in Supplementary File 1.

Understanding gene expression, as already described, requires analysis of kinetic and structural information pertaining to gene regulatory interactions. In the absence of the former, we model gene expression as a signal propagation process [START_REF] Anchang | Modeling the temporal interplay of molecular signaling and gene expression by using dynamic nested effects models[END_REF][START_REF] Stetter | Large-scale computational modeling of genetic regulatory networks[END_REF], by which a regulatory gene attempts to modify the activity of particular target genes. Therefore, we model the signal propagation, itself, as a random walk on a directed graph. Diffusion processes (such as propagation) have already been considered and modeled via diffusion kernels, mostly with focus on undirected graphs [START_REF] Kondor | Diffusion kernels on graphs and other discreet structures[END_REF][START_REF] Simonsen | Diffusion on complex networks : A way to probe their large scale topological structures[END_REF].

For a vector x, let diag(x) denote the diagonal matrix, which has the entries of vector x on its diagonal, and all other entries equal zero. Given a directed graph G = (V, E), a random walk on G with adjacency matrix A is a Markov process with transition matrix P = D -1 A, where D = diag(Ae) is a diagonal matrix with entries on the diagonal corresponding to the out-degrees of the nodes in V (G) (e is the vector in which all entries are 1). Assume, for now, that G is strongly connected i.e., for every u, v ∈ V (G), there exist directed paths, from u and v and from v to u. The Perron-Frobenius theorem for matrix P then states that there exists a unique left eigenvector which is strictly positive with eigenvalue 1 (since P e = e) [START_REF] Ninio | A simple proof of the Perron-Frobenius theorem for positive symmetric matrices[END_REF][START_REF] Seneta | Non-negative matrices and Markov chains[END_REF]. Note that P has a unique normalized left eigenvector with eigenvalue 1 if G is aperiodic [START_REF] Langville | Deeper inside PageRank[END_REF].

Let π be the unique normalized left eigenvector such that

πP = π (1)
and u∈V (G) π(u) = 1. Furthermore, the row-vector π corresponds to the stationary distribution of the random walk defined by P . From the Eq. ( 1), we have that

π(u) = v,v→u π(v)P (v, u),
(2)

i.e., the probability of finding the random walk at u is the sum of all incoming probabilities to nodes v which are neighbours of u. We can now define a circulation of the directed graph G.

Definition 0.1. [START_REF] Chung | Laplacians and the cheeger inequality for directed graphs[END_REF] A function

F : E(G) → R + 0 that assigns each directed edge to a non-negative value is called a circulation if u,u→v F (u, v) = w,v→w F (v, w),
for each node v.

The circulation can be easily interpreted as a flow in the graph. The flow at each node must be conserved, hence, the flow in is equal to the flow out. One such circulation is defined in terms of the stationary distribution of the random walk on G (see Eq.

(2)) and is given by

F π (u, v) = π(u)P (u, v).
If the random walk is reversible, then

F π (u, v) = F π (v, u).
We can now examine the matrix:

à = ΠP + P T Π 2 , ( 3
)
where Π = diag(π). In contrast to A, the adjacency matrix of the directed graph G, the matrix Ã, is symmetric and, thus, has an undirected graph G corresponding to G. In contrast to G, which is unweighted, the obtained undirected graph G is edge-weighted.

The weights of the edges in G correspond to the defined circulation in terms of the random walk on G.

The Laplacian of the graph G is L( G) = Π -Ã. To eliminate the influence of the number of nodes and edges, we investigate the normalized directed graph Laplacian, which is defined in terms of L( G) as follows:

L(G) = Π -1/2 Π -Ã Π 1/2 = I - 1 2 Π 1/2 P Π -1/2 + Π -1/2 P T Π 1/2 . ( 4
)
Initially, we assumed that the graph G is strongly connected and aperiodic, and therefore that G has a unique stationary distribution. If G is not strongly connected, one can define the Page rank transformation matrix, which is a transformation of the transition matrix P to ensure aperiodicity and strong connectivity (i.e., there is a path between any two nodes) [START_REF] Bianchini | Inside PageRank[END_REF][START_REF] Langville | Deeper inside PageRank[END_REF]. The transformation can be formalized as follows:

P pr = αP + (1 -α) n ee T , ( 5
)
where α = 0.85, to ensure convergence of the random walk. An example of the transformation imposed by Eq. ( 5) appears in Supplementary File 2.

Although the modifications by Eqs.

(3) and ( 5) alter the directed graph G, when applied in the same fashion across different directed graphs, they allow for systematic network comparison. One way to compare different directed graphs is through the spectra of their Laplacians. The spectrum of L(G) is given by n eigenvalues

λ 0 = 0 ≤ λ 1 ≤ . . . ≤ λ n-1 .
If the graph G is strongly connected, then λ 1 > λ 0 = 0 and λ 1 is known as the algebraic connectivity of the directed graph G.

The algebraic connectivity λ 1 of a graph G is connected to another graph invariant called Cheeger constant [START_REF] Chung | Spectral Graph Theory[END_REF]. The Cheeger constant is a numerical measure of whether or not a graph has a "bottleneck". To define the Cheeger constant for a directed graph, we need to introduce the concept of out-boundary: Definition 0.2. Given a graph G and a subset of nodes S ⊂ V (G), the out-boundary of S, denoted by ∂S, is comprised of all edges going from a node in S to a node in S,

S = V (G) -S, i.e., ∂S = {(u, v) ∈ E(G)|u ∈ S, v ∈ S}.
Definition 0.3. The volume of a node u is defined as vol(u) = v P (u, v), the volume of set S is defined as vol(S) = u∈S vol(u), and the volume of an out-boundary is defined as vol(∂S) = u∈S,v∈ S P (u, v).

The Cheeger constant for a graph G, denoted by h(G), can then be defined in terms of the Cheeger constant of a cut (S, S) implied by a given node-subset S ⊂ V (G):

Definition 0.4. The Cheeger constant of a cut (S, S) in a graph G is h G (S) = vol(∂S) min{vol(S), vol( S)} .
The Cheeger constant of the graph G is

h(G) = min S⊂V (G) h G (S).
The Cheeger constant is strictly positive if and only if G is a connected graph. Intuitively, if the Cheeger constant is small but positive, then there exists a "bottleneck", in the sense that there are two large sets of nodes with few edges between them. The Cheeger constant is large if any possible bi-partition of the node set has many edges between the two partitions. Since computing the Cheeger constant is NP-hard, we rely on its approximation in terms of the algebraic connectivity expressed by the following theorem:

Theorem 0.5. [START_REF] Chung | Spectral Graph Theory[END_REF] The algebraic connectivity, λ 1 , of the normalized directed

Laplacian of the graph G is related to the Cheeger constant h(G) by

h(G) 2 2 ≤ λ 1 ≤ 2h(G).
With these notions, we can now analyze the spectrum of the directed TF-gene and gene-gene regulatory networks via the generalization of graph Laplacian for directed graphs. As explained, this generalization symmetrizes the graph (i.e., turns it into an undirected graph) without loss of information. From the directed graph Laplacian, we determine the algebraic connectivity, and establish connections between this graph invariant for both, TF-gene and gene-gene, regulatory networks and their corresponding average path length (from the distribution of cascades), diameter, and average degree. In addition, we determine the local structures shown to have effect on the dynamics of gene expression, namely the concentrations of all 2-and 3-node directed subgraphs, and then calculate their Pearson correlation with the algebraic connectivity across the analyzed organisms. The concentration of a subgraph i on n nodes occurring M n i was calculated by

c n i = M n i i M n i .
Based on the correlation we can rank the importance of a particular property and local structure in the evolution of regulatory architectures. Note that we do not assess the significance of a subnetwork, since we would like to determine how the abundance of each of these subnetworks correlates with the algebraic connectivity.

In addition, we calculate the evolutionary tree from the considered network properties by first determining the distance matrix and then using it to perform hierarchical clustering. We define the distance between two organisms with respect to a network property as the absolute difference between the values of the property for the networks of the two organisms. This is the simplest distance measure which already shows results in good agreement with the established phylogenetic relationships of organisms analyzed here.

For distance-based methods for building phylogenetic trees from metabolic networks the reader is directed to [START_REF] Mazurie | Phylogenetic distances are encoded in networks of interacting pathways[END_REF] and reference therein.

Data set

We analyze the data set comprising the interactions between regulatory genes, transcription factors, and target genes, extracted from TRANSFAC 7.0 database [START_REF] Matys | Transfac: transcriptional regulation, from patterns to profiles[END_REF], for 97 organisms from all kingdoms of life. TRANSFAC contains data on transcription factors, their experimentally proven binding sites, and regulated genes. In the present study, an edge between a gene and a transcription factor is established if the gene encodes the transcription factor; moreover, an edge is established between a transcription factor and a gene if the transcription factor is experimentally proven to bind to the gene.

In addition to the networks from TRANSFAC we analyzed curated gene regulatory networks for Escherichia coli, extracted from Regulon DB and the developmental gene regulatory networks of Sea urchin's ectoderm and endomesoderm, obtained from [START_REF] Davidson | A genomic regulatory network for development[END_REF]. In total, our preliminary analysis is conducted on 100 gene regulatory networks. The analysis was performed by using the igraph R package.

The assembled TF-gene and gene-gene regulatory networks were first subjected to preprocessing in order to identify the weakly connected components (WCCs) and their corresponding distributions across all organism-specific networks. From the 100 analyzed TF-gene regulatory networks, 17 have a largest weakly connected component which contains more than 40% of the nodes from the originally assembled network. The other connected components do not contain more than 4 nodes, and were excluded from further consideration to reduce the bias due to different number of small WCCs. From the 17 identified networks, only 9 contain more than 80 nodes (see Supplementary File 3). The remaining 8 networks contain small number of nodes, between 4 and 20, and were also excluded from further analysis. Finally, our analysis was performed on the largest WCCs of 8 networks from the following 7 organisms: Escherichia coli, Saccharomyces cerevisiae, Sea urchin (two networks-ectoderm and endomesoderm), Drosophila melanogaster, Mus musculus, Rattus norvegicus, and Homo sapiens.

In the following, we present the results for the properties of the largest weakly connected component of the 8 organism-specific TF-gene and gene-gene regulatory networks.

The biological interpretations of the obtained results ultimately depend on the quality and scope of the considered data. Therefore, the application of the proposed method on these and other data sets warrants revisiting the analysis in a future study when the quality of known gene-regulatory networks may be considerably improved.

Results and Discussion

Network properties of TF-gene and gene-gene regulatory networks

The average degree of the TF-gene regulatory networks was found to be in the range from 1.235 for Saccharomyces cerevisiae to 4.821 for Sea urchin ectoderm (see Table 1).

Since for any directed graph, the average in-and out-degrees are the same, we report only one value. There is an apparent difference in the average degree between TF-gene networks over the entire genome and regulatory networks involved in developmental pro- to gene-gene networks. For the gene-gene regulatory networks, the average degree ranges from 1.322 for Saccharomyces cerevisiae to 6.836 for Sea urchin ectoderm, as shown in Table 2. However, like for the TF-gene networks, the evolutionary relationships cannot be recovered from this network property, as shown in Fig. 1 (A).

Previous analyses have demonstrated differences in the distribution of cascades between Escherichia coli and Drosophila melanogaster, and speculated the biological implications of this finding [START_REF] Rosenfeld | Response delays and the structure of transcription networks[END_REF]. Namely, [START_REF] Rosenfeld | Response delays and the structure of transcription networks[END_REF] found that the longest cascade in the network of Escherichia coli is of length 4, and that the majority of cascades are of unit length; in the network of Drosophila melanogaster, the largest cascade is of length 9 and the majority of cascades are of length between 6 and 9. Our results for the distribution of cascade lengths in TF-gene regulatory networks, analyzed here, are shown in Fig. 2 for Escherichia coli closely follows those of Sea urchin's networks. These findings point out that the distribution of cascade lengths may not play a significant role in shaping the differences between the investigated organisms.

The results for gene-gene regulatory networks differ by a factor of 2, due to the transformation from TF-gene to gene-gene networks, which turns every directed path of length 2 into a directed edge, as shown in Table 2. From Fig. 1 (B) and (C), we can conclude that neither the length of the longest cascade nor the average cascade length can be used to recover the evolutionary relationship between the investigated organisms.

A typical result of discordance is the clustering of Escherichia coli and Homo sapiens, in the case of the longest cascade, and Homo sapiens and Sea urchin, in the case of the average path length.

Connected subnetworks in gene-gene regulatory networks

The distribution of connected subnetworks on n nodes can be regarded as a feature of a gene-gene regulatory network. As the number of nonisomorphic directed connected graphs on n nodes is a fast growing function [START_REF] Harary | Graphical Enumeration[END_REF], we determined the concentrations for each of 2-and 3-connected subnetworks. As in [START_REF] Milo | Network motifs: Simple building blocks of complex networks[END_REF], a subnetwork was considered to be significant if it appears at least 4 times in an analyzed network. Note that we are not interested in establishing the significance of the concentration of a particular subnetwork, but only on the possible correlation between the subgraph concentration and other structural properties. The results for the distribution of 3-connected subnetworks appear in Fig. 3.

The two 2-connected subnetworks, i.e., a directed edge and a pair of directed edges in the investigated networks, cannot discriminate the complexity of the organisms and behave similarly to the average degree. The significant subnetworks across all organisms include the following: M1, M2, M4, and M5 (see Fig. 3). The subnetworks M3, M4, M6, and M11 are significant for all organisms except Saccharomyces cerevisiae and Rattus norvegicus, while subnetworks M8, M9, M10, M12, and M13 are not significant in all networks except the two networks of Sea urchin. While the last five enumerated subnetworks do separate the developmental regulatory networks from those on the entire genome, we still face the same problem as in the case of the distribution of cascade lengths:

Namely, the distributions of subnetwork concentrations of Saccharomyces cerevisiae and Homo sapiens are closest to each other and the distribution of Drosophila melanogaster is closest to those of Sea urchin's networks.

Algebraic connectivity of gene-gene regulatory networks

To account for the various endogenous and exogenous conditions which may affect the expression of genes in a given organism, we model gene expression as signal propagation,

i.e., a random walk by which a regulatory gene can affect a target gene. The convergence of this random walk to its stationary distribution is characterized by the algebraic connectivity of the directed Laplacian of the graph on which the random walk takes place.

The algebraic connectivity is related to the Cheeger constant (see Methods) which characterizes the existence of bottlenecks in the underlying network structure. The bottlenecks can be seen as the means for throttling or slowing down signal propagation on the graph.

In graphs with large Cheeger constant, signals can easily propagate between any two nodes, since any two node partitions are connected with a large number of edges. Such a property is desirable for gene regulatory networks which must settle in a stationary state in a shorter time, rendering, ultimately, shorter response times on a system's level. On the other hand, if a subset of regulatory interactions is condition-specific, its connection to the rest of the network should be limited, and only a portion of the network is affected (i.e., the random walk, on average, remains localized in one subnetwork).

The seemingly absent link between algebraic connectivity of TF-gene networks and evolutionary relatedness, which follows from Table 1, is due to the fact that TF-gene networks are bipartite . As stated in the seminal paper of [START_REF] Lovasz | Random walks on graphs[END_REF], the Markov chain P does not converge to a stationary distribution. Moreover, the PageRank algorithm used in the derivation of the algebraic connectivity has the same degeneracy on bipartite graphs, yielding different ranking for nodes of equal indegree (see Theorem 4 [START_REF] Meghabghab | Search Engines, Link Analysis, and Users Web Behavior: A Unifying Web Mining Approach[END_REF]). This drawback of using TF-gene networks is overcome precisely by the transformation of the bipartite TF-gene graphs into non-bipartite gene-gene graphs.

The algebraic connectivity of the gene-gene regulatory networks demonstrated an in-teresting property-the more complex the organism, the smaller the algebraic connectivity (see Table 2). Therefore, organisms whose gene regulatory interactions should react with higher specificity to different stimuli depending on the tissue/conditions, exhibit smaller algebraic connectivity, effectively rendering localization of the random walk (e.g., the networks of Drosophila melanogaster, Rattus norvegicus, Mus musculus, and Homo sapiens have algebraic connectivity of 0.009, 0.007, 0.009, and 0.007, respectively). Moreover, the developmental networks of Sea urchin have high algebraic connectivity (0.109 for ectoderm and 0.010 for endomesoderm), comparable to that of Escherichia coli (0.08) and

Saccharomyces cerevisiae (0.02). This implies that the gene regulatory interactions are structured to facilitate faster convergence of gene expression to a stationary state for organisms and processes that should adapt faster to exogenous conditions. Furthermore, as shown in Fig. 1 (D), the value of the algebraic connectivity does reflect the evolutionary relationship between the investigated organisms.

The Pearson correlation between the average degree and the algebraic connectivity over the investigated organisms is 0.640, which implies that the average degree does have influence on the algebraic connectivity. Recent mathematical result points at the rules for building a graph model for which the algebraic connectivity can be made arbitrarily close regardless of the given (prescribed) degree distribution [START_REF] Atay | Synchronization of networks with pre-scribed degree distributions[END_REF]. Therefore, our finding suggests that, although the degree distribution of the regulatory networks from different organisms may be scale-free (even with equal exponents, Fig. 4 (A) and (B)), the edges in the graph are distributed in a way that the increased complexity is matched with smaller algebraic connectivity.

On the other hand, the correlation between the longest cascade and the algebraic connectivity over the investigated organisms is -0.19, which shows that the longest cascade has only a negligible reciprocal effect on the algebraic connectivity. The same can be observed for the correlation between the average path length and algebraic connectivity.

Recent theoretical results demonstrated that synchronizability of the complex networks worsens with the increase of the average path length [START_REF] Zhao | Relations between average distance, heterogeneity and network synchronizability[END_REF]. Since the Laplacian matrix arises naturally in the study of models of synchronization [START_REF] Watts | Six degrees: The science of a connected age[END_REF], our empirical results point out that gene regulatory networks form a class of complex networks whose synchronizability may not be affected by the average path length, but is tightly coupled with the algebraic connectivity.

Next, we determine to what extent the local structures determine the existence of the bottlenecks which may have effects on gene expression. To this end, we calculated the correlation between the 3-connected subnetworks (see Fig. 3) and algebraic connectivity across the investigated organisms. The algebraic connectivity exhibits the largest correlation of 0.934 with the concentration of subnetworks M6, followed by M5 (0.886), M8 (0.769), M10 (0.766), M9 and M12 (0.766). From these ranking of subnetworks and the previous results about their significance, we can conclude that subnetworks M5 and M6

have the highest influence on the value of the algebraic connectivity.

Although, in gene-gene networks, the concentrations of subgraphs M5 and M6, see Table 2, show hight correlation with the algebraic connectivity, only the latter reflects the evolutionary relatedness between the considered organisms (see Supplementary File 4). The reason for the discrepancy in the evolutionary relatedness according to the two properties is due to the differences in the type of structural information employed for calculating them. For instance, the concentration of a particular 3-subnetwork is calculated based on a local property of how 3 nodes are connected to each other. Therefore, it can only capture a limited portion of the network topology. Moreover, the concentration of a subgraph states neither how the copies of the subnetwork are distributed across the network, nor how they are positioned in conjunction to the remaining types of sunet-works. Hence, deriving principles of network evolution on a local property may represent an ill-founded approach, as already illsutrated.

On the other hand, algebraic connectivity is a global property, whose calculation depends on global information captured in the matrices involved in its calculation. It is likely that the algebraic connectivity, as a global property, reflects the information about not only the concentration of a particular subnetwork, but also of how the copies of the subnetwork are distributed and positioned in the network. Making a formal statement and empirical investigation about the relation of the algebraic connectivity and the concentration of subnetworks in complex networks is far beyond the scope of this paper, and remains as an interesting open problem.

A closer look at the subnetworks M5 and M6, identified to have the highest rank with respect to their correlation with the algebraic connectivity, we conclude that M5 is precisely the feed forward loop (FFL) and M6 contains the FFL as an embedded subnetwork.

The feed forward loop is determined by three genes, say u, v, and w, such that u regulates v, while both u and v regulate the gene w. The FFL has been demonstrated theoretically and experimentally to perform a basic information-processing function, i.e., it shows a delay following ON steps of an input inducer, but not after OFF steps. Moreover, it was recently demonstrated that the FFL is selected over simpler subnetworks in environments where the distribution of the input pulse duration is sufficiently broad and contains both long and short pulses [START_REF] Dekel | Environmental selection of the feed-forward loop circuit in gene-regulation networks[END_REF][START_REF] Mangan | Structure and function of the feed-forward loop network motif[END_REF]. In addition, recent simulation studies have demonstrated that in directed graphs, synchronization (directly expressed via the algebraic connectivity) can be correlated with the abundance of FFL as a specific 3-node subnetworks [START_REF] Brede | Synchronization on directed small wolrds: Feed forward loops and cycles[END_REF][START_REF] Lodato | Synchronization properties of network motifs[END_REF].

With regard to the ranking of subnetworks, our result on the positive correlation between the concentration of subnetworks M5 and the algebraic connectivity implies that the smaller algebraic connectivity requires smaller concentration of strategically embedded subnetworks isomorphic with the FFL. Therefore, in addition to the identified function of single FFL to detect non-temporary change of the environment, our results suggest that evolution favours decreasing concentration of FFLs combined in such a way that the network, as a whole, can throttle a random walk on the network and, thus, increase the specificity towards changing environments.

Our results point out that the algebraic connectivity is a salient property of genegene regulatory networks and encodes other network properties (e.g., average degree and subnetwork concentration) governing the network evolution. Although the FFL subnetwork and scale-free degree distributions have already been identified as salient properties of gene regulatory networks, it was not known what is the evolutionary principle which select these over other network properties. Moreover, thus far, it was not known why the FFLs appear as significant in virtually all gene regulatory networks as well as their relations with regulatory networks coding for different organismal complexity. Here, we identified the decrease in the algebraic connectivity of the directed regulatory network as one possible evolutionary principle.

Although at present we cannot exhibit a model to mimic these findings, our results point out at two necessary rules for developing gene regulatory networks which encode growing complexity: (1) the degree distribution has a heavy-tail with average degree positively correlated with the decreasing algebraic connectivity, and (2) the distribution of edges is such that the concentration of the FFL subnetwork is positively correlated with the decreasing algebraic connectivity. Initial steps for developing the first rule on undirected graphs have already been made in [START_REF] Atay | Synchronization of networks with pre-scribed degree distributions[END_REF]. The second rule can be addressed by developing a model of directed networks which considers the effect of gene duplication events and subsequent sub-and neo-functionalization of genes on the

  cesses. If the average degree in both types of networks played a role in the evolution of these regulatory networks, one would not have expected the similar behavior of this network property on networks from organisms of different kingdoms, such as: Saccharomyces cerevisiae (1.235) and Mus musculus (1.370) or Escherichia coli (2.122) and Sea urchin endomesoderm (3.588) with closest values to each other. The difference in the average degree across organisms becomes more pronounced for the gene-gene regulatory networks due to the densification implied by the transformation from bipartite TF-gene

Figure 1 .

 1 Figure 1. Evolutionary relationship from network properties of gene-gene networks. Heatmaps together with the corresponding hierarchical clustering based on: (A) average degree, (B) longest cascade, (C) average cascade length, and (D) algebraic connectivity for Escherichia coli (Ec), Saccharomyces cerevisiae (Sc), Sea urchin ectoderm (SUect), Sea urchin endomesoderm (SUend), Drosophila melanogaster (Dm), Mus musculus (Mm), Rattus norvegicus (Rn), and Homo sapiens (Hs).

Figure 2 .

 2 Figure 2. Distribution of cascade lengths in TF-gene networks across organisms. The distribution for Saccharomyces cerevisiae (green line) confirms its shallow architecture, with Rattus norvegicus being its closest neighbour. Moreover, the distribution of cascades between Escherichia coli and Drosophila melanogaster differ only slightly, in contradiction to already reported results.

Figure 3 .

 3 Figure 3. Distribution of 3-connected subnetworks in gene-gene networks of studied organisms. The thirteen 3-connected directed subnetworks are shown on the right. Subnetwork M5 is the feed forward loop. Two groups of organisms can be observed based on the correlation between the distributions of 3-connected subnetworks: the first is given by Drosophila melanogaster and Sea urchin, the second by Homo sapiens, Mus musculus, Rattus norvegicus, and Saccharomyces cerevisiae. The network of Escherichia coli stands on its own, due to the large concentration of the subnetwork M1.

Figure 4 .

 4 Figure 4. Degree distributions. Distribution of the (A) in-degree and (B) out-degree of the analyzed networks. All networks, except for those of Sea urchin, follow the same trend in the degree distribution for both in-and out-degrees.

  . The shape of the distributions for Saccharomyces cerevisiae and Rattus norvegicus point out at their shallow regulatory architectures, with predominantly small length cascades. The longest cascades are found in Homo sapiens (17), Sea urchin endomesoderm (14), and Escherichia coli (12). Moreover, the shape of the distribution of cascade lengths in Escherichia coli (red) and Drosophila melanogaster (blue) are very similar, as shown in Fig.2, which is contrary to the previously reported

results. Furthermore, the average path length is almost identical for Escherichia coli (3.804) and Drosophila melanogaster (3.884). The largest average path lengths are found in Homo sapiens (5.249) and Sea urchin endomesoderm (3.973), followed by Drosophila melanogaster and Escherichia coli. Finally, the shape of the cascade length distribution

Table 1 . Network properties for TF-gene networks

 1

	Tables

Network properties of TF-gene regulatory networks.

  Six network properties: number of nodes (n), number of directed edges (m), algebraic connectivity (λ 1 ), longest cascade (D), average cascade (D avg ), and average degree (d avg ) for TF-gene regulatory network of Escherichia coli, Saccharomyces cerevisiae, Sea urchin's ectoderm and endomesoderm, Drosophila melanogaster, Mus musculus, Rattus norvegicus, and Homo sapiens.

Table 2 . Network properties of gene-gene networks

 2

	network properties

Network properties of gene-gene regulatory networks.

  Six network properties: number of nodes (n), number of directed edges (m), algebraic connectivity (λ 1 ), longest cascade (D), average cascade (D avg ), average degree (d avg ), concentration of M5 subnetwork (c M 5 ), and concentration of M6 subnetwork (c M 6 ) for gene-gene regulatory network of Escherichia coli, Saccharomyces cerevisiae, Sea urchin's ectoderm and endomesoderm, Drosophila melanogaster, Mus musculus, Rattus norvegicus, and Homo sapiens.

Acknowledgments

Z.N., P.M., and J.S. would like to acknowledge the support from the GoFORSYS project funded by the German Federal Ministry of Education and Research, Grant Nr. 0313924.

algebraic connectivity.