text:
Human blood myeloid and plasmacytoid dendritic cells cross activate each other and synergize in inducing NK cell cytotoxicity

ABSTRACT Human blood dendritic cells (DCs) hold great potential for use in anticancer immunotherapies. CD1c+ myeloid DCs and plasmacytoid DCs (pDCs) have been successfully utilized in clinical vaccination trials against melanoma. We hypothesize that combining both DC subsets in a single vaccine can further improve vaccine efficacy. Here, we have determined the potential synergy between the two subsets in vitro on the level of maturation, cytokine expression, and effector cell induction. Toll-like receptor (TLR) stimulation of CD1c+ DCs induced cross-activation of immature pDCs and vice versa. When both subsets were stimulated together using TLR agonists, CD86 expression on pDCs was increased and higher levels of interferon (IFN)-α were produced by DC co-cultures. Although the two subsets did not display any synergistic effect on naive CD4+ and CD8+ T cell polarization, CD1c+ DCs and pDCs were able to complement each other's induction of other immune effector cells. The mere presence of pDCs in DC co-cultures promoted plasma cell differentiation from activated autologous B cells. Similarly, CD1c+ DCs, alone or in co-cultures, induced high levels of IFN-γ from allogeneic peripheral blood lymphocytes or activated autologous natural killer (NK) cells. Both CD1c+ DCs and pDCs could enhance NK cell cytotoxicity, and interestingly DC co-cultures further enhanced NK cell-mediated killing of an NK-resistant tumor cell line. These results indicate that co-application of human blood DC subsets could render DC-based anticancer vaccines more efficacious.

Introduction
Dendritic cells (DCs) possess the ability to drive both adaptive and innate immune responses and therefore play a central role in antitumor immunity. Hence, they can be exploited for the treatment of immunogenic cancers, such as melanoma. 1,2 DC-based vaccines consist of patient-derived DCs, which are stimulated and loaded with tumor-associated antigens (TAAs) ex vivo, before being reinfused into the patient. This form of immunotherapy aims at activating the patient's own immune system to recognize and eradicate tumor cells. DCs derived from monocytes or CD34 C progenitor cells are usually the main source of DCs used for vaccines.
Nevertheless, the use of naturally occurring DCs isolated from blood has been recently hypothesized to improve therapeutic responses, due to less extensive ex vivo culture periods. 3 Several DC subsets can be identified in human peripheral blood. 4 They are divided into plasmacytoid DCs (pDCs) and myeloid DCs (mDCs), with the latter being subdivided into CD1c C (or BDCA1 C ) DCs and the rare CD141 C (or BDCA3 C ) DCs. mDCs and pDCs are functionally distinct, which is reflected by their non-overlapping repertoire of expressed tolllike receptors (TLRs). [5][6][7] Cytokine secretion following stimulation also differs between DC subsets. CD1c C DCs can secrete high levels of bioactive interleukin (IL)-12, 8 an important cytokine for the induction of T helper 1 (Th1) and cytotoxic T lymphocyte (CTL) responses. 9,10 In contrast, pDCs can produce massive amounts of type I interferons (IFNs) upon stimulation. 11 IFN-a can also take part in the skewing of Th1 responses, 12 in addition to increasing the cytotoxic activity of natural killer (NK) cells. 13 These characteristics make both DC subsets suitable for use in cancer immunotherapy. Indeed, vaccination with CD1c C DCs for prostate cancer was shown to be safe and feasible. 14 Our group has previously conducted phase I clinical trials exploiting either pDCs or CD1c C DCs for vaccination of melanoma patients, demonstrating the safety and efficacy of this approach. 15,16 These promising results raise the question whether combining both DC subsets can further improve immunological responses and clinical outcome.
In the current study, we have characterized the crosstalk between human blood CD1c C DCs and pDCs by analyzing maturation status and cytokine production after co-culture of both subsets in the presence or absence of TLR stimulation. In addition, functional implications of crosstalk-mediated effects on other adaptive and innate immune cells, namely naive CD4 C and CD8 C T cells, B cells and NK cells, were addressed. We have shown that CD1c C DCs and pDCs are able to cross activate each other. Although co-application of DC subsets did not augment T cell polarization, it did allow for complementation of subset-specific effector functions, including induction of plasma cell differentiation from B cells and secretion of high levels of IFN-g by peripheral blood lymphocytes (PBLs) and NK cells. Furthermore, DC co-cultures could synergistically enhance NK cell-mediated killing of an NK-resistant tumor cell line. These results indicate that combining the distinct qualities of different human blood DC subsets may potentiate current DC-based anticancer vaccines for optimal therapeutic outcome.

TLR triggering of human CD1c C DCs or pDCs cross activates bystander DCs
We first investigated whether human DC subsets can directly alter each other's activation status. Autologous human CD1c C DCs and pDCs were isolated from peripheral blood and stimulated overnight, either separately (single cultures) or together (co-cultures) at a 1:1 ratio. The total number of DCs was kept constant between the cultures. DCs were stimulated with a specific TLR ligand: CpG oligodeoxynucleotide (ODN) class C (CpG-C) for pDCs, polyinosinic-polycytidylic acid (poly(I:C)) for CD1c C DCs. In the absence of TLR agonists, IL-3 and granulocyte-macrophage colonystimulating factor (GM-CSF) were added to maintain the survival of pDCs and CD1c C DCs, respectively. Subset-specific maturation status was determined by flow cytometry analysis of the co-stimulatory molecules CD86, CD80, and CD40, the maturation marker CD83, the co-inhibitory molecule programmed deathligand 1 (PD-L1), C-C chemokine receptor type 7 (CCR7), and major histocompatibility complex (MHC) class I and II molecules. The presence of poly(I:C)-activated CD1c C DCs led to significantly enhanced expression levels of CD86 (p D 0.0243), PD-L1 (p D 0.0002) and MHC class I (p < 0.0001) on unstimulated pDCs (Fig. 1A). Conversely, immature CD1c C DCs significantly upregulated CD86 (p D 0.0445), CD80 (p < 0.0001), PD-L1 (p < 0.0001), and MHC class I (p D 0.0014) and slightly downregulated MHC class II (p D 0.0414), when co-cultured with CpG-Cactivated pDCs (Fig. 1B). When both subsets were stimulated together, using their respective TLR agonists, pDCs showed increased CD86 expression compared to CpG-C-stimulated pDCs cultured separately (p D 0.0464), whereas CD1c C DCs showed a trend toward increased PD-L1 expression (p D 0.0741). To exclude the possibility that changes in DC maturation status were the result of a combination of stimuli, rather than crosstalk between DC subsets, combinations of stimuli that were used for DC co-cultures were also tested on DC single cultures. Incorporation of these extra control groups did not change the above mentioned results qualitatively ( Fig. S1 A and B). Of note, the combination of CpG-C and poly(I:C) did not fully mature CD1c C DCs in both DC single cultures and co-cultures, as co-stimulatory molecules were upregulated less when poly(I:C) was used as the sole stimulus.
Next, cytokine secretion by human blood DCs was examined. Poly(I:C)-stimulated CD1c C DCs secreted IL-12p70, the bioactive form of IL-12, although no detectable levels were reached when a combination of poly(I:C) and CpG-C was used (Fig. 1C, first panel and Fig. S1C). In contrast, CpG-C-stimulated pDCs produced high levels of IFN-a. Interestingly, IFN-a secretion was significantly augmented in DC co-cultures where CpG-C and poly(I:C) were used, compared to CpG-Cstimulated pDCs alone (p D 0.0255; Fig. 1C, second panel). Furthermore, both subsets produced tumor necrosis factor (TNF)-a upon TLR ligation. When poly(I:C)-stimulated CD1c C DCs were co-cultured with immature pDCs, there was a trend toward more TNF-a secretion compared to the corresponding DC single cultures (p D 0.1537 when compared to poly(I:C)treated CD1c C DCs alone; Fig. 1C, third panel). These results indicate that TLR triggering of one DC subset induces cross activation of the other subset in culture. In addition, TLR triggering of both DC subsets may induce further changes to their phenotype and functionality, such as increased expression of co-stimulatory molecules and secretion of cytokines.
To assess whether paracrine cytokine stimulation is the main mechanism behind the upregulation of maturation markers observed in the DC co-cultures, immature CD1c C DCs were cultured with different concentrations of recombinant IFN-a or TNF-a, whereas immature pDCs were cultured with different combinations of recombinant TNF-a. IFN-a by itself induced strong upregulation of CD86, CD80, PD-L1, and MHC class I on CD1c C DCs (Fig. S2A), all to a similar extent as can be seen in the DC co-cultures (Fig. 1). Incubation with TNF-a alone induced upregulation of PD-L1 and MHC class I and a slight decrease of CD86 on immature CD1c C DCs (Fig. S2B), whereas at the highest tested concentration, it induced a slight but significant upregulation of CD80 on immature pDCs (p D 0.0253; Fig. S2C). Thus, although TNF-a plays a partial role in the crosstalk between pDCs and CD1c C DCs, the cross activation of CD1c C by activated pDCs seems predominantly mediated by IFN-a.
Combined application of CD1c C DCs and pDCs sustains their capacity to activate lymphocytes Following our observation that CD1c C DCs and pDCs could cross activate each other, we investigated whether this may influence the DCs' ability to activate various immune effector cells. Since poly (I:C) did not optimally mature CD1c C DCs in the presence of CpG-C, stimulation by lipopolysaccharide (LPS) was included as a condition in all further assays. LPS induced strong maturation of mDCs, but not pDCs, while no significant IL-12 production was induced (data not shown). First, a mixed lymphocyte reaction (MLR) assay was performed to determine if DC co-cultures have an effect on the proliferation and cytokine production of PBLs. Expectedly, both CD1c C DCs and pDCs induced proliferation of allogeneic PBLs. However, the combination of DC subsets did not have an additive effect on PBL proliferation (Fig. 2, first panel). Furthermore, both DC subsets induced IFN-g secretion in the MLR, although higher IFN-g levels were reached with TLR-stimulated CD1c C DCs than with pDCs ( Fig. 2, second panel). Similar to PBL proliferation, DC co-cultures did not yield superior IFN-g production by PBLs. Collectively, these results clearly demonstrate that human DC subsets, when cultured together, maintain their capacity to induce PBL activation, yet without displaying any synergy.
Co-cultures of CD1c C DCs and pDCs drive CD4 C and CD8 C T cell polarization DCs are important for the polarization of CD4 C T cells into the different T helper cell subsets, as well as the induction of CTLs from naive CD8 C T cells. Tumor-specific CD8 C T cell responses are critical for the eradication of tumor cells in cancer immunotherapy. In addition, tumor-specific CD4 C T cells contribute to CD8 C T cell responses and their targeting enhances efficacy of DC vaccination for treatment of melanoma patients. 17,18 For this purpose, these effector CD4 C T cells should exhibit a Th1 phenotype, which is characterized by secretion of IFN-g. 19 In sharp contrast, differentiation of CD4 C T cells into IL-10-secreting regulatory T cells inhibits adaptive immune responses and thereby promotes tumor immune evasion. 20 To determine if a combination of human blood DC subsets can efficiently polarize T cells into effector cells that are favorable for cancer immunotherapy, CD1c C DCs and/or pDCs were used to activate allogeneic naive CD4 C or CD8 C T cells. T cell polarization was determined by analyzing cytokine production profile of the induced T cells. Both DC single cultures and DC co-cultures induced IFN-g C CD4 C T cells (Fig. 3A, first panel). Total IFN-g secretion in the CD4 C T cell cultures differed slightly between the different DC stimuli, but no significant changes were detected for the combination of CD1c C DCs and pDC, compared to the use of either subset alone (Fig. 3A, second panel). In addition, a favorable ratio of IFN-g to IL-10 was induced in DC co-cultures and IL-10 production in DC co-cultures was lower than in cultures with CpG-C-activated pDCs (Fig. 3A, third panel). These results indicate that the combination of CD1c C DCs and pDCs predominantly drives Th1 polarization, with low induction of regulatory T cells. Finally, CD1c C DCs, pDCs or the combination of both induced granzyme B C CD8 C T cells, IFN-g C CD8 C T cells and TNF-a C CD8 C T cells, with no significant advantage of DC cocultures over DC single cultures (Fig. 3B). Thus, combining CD1c C DCs and pDCs does not further augment T cell polarization in comparison to the application of either DC subset alone.

Co-cultures of CD1c C DCs and pDCs promote plasma cell differentiation of activated B cells
Although combined application of pDCs and CD1c C DCs did not have an added value in T cell activation, we investigated whether it has an effect on other immunological components that are essential for tumor eradication. B cells have a crucial role in antitumor immunity as antibody-secreting plasma cells. 21 Plasma cell differentiation from naive B cells can be promoted by activated pDCs. 22 We examined if this effect persists in co-cultures of human CD1c C DCs and pDCs. Freshly isolated CD1c C DCs and/or pDCs were incubated together with autologous B cells in the presence of CpG-C and other stimuli for 7 d. CpG-C acts as a stimulus for both pDCs and B cells, and poly(I:C) or LPS was added to stimulate CD1c C DCs. In sharp contrast to CD1c C DCs, activated pDCs enhanced plasma cell differentiation of activated B cells ( Fig. 4 A and B). Combining DC subsets did not enhance the plasma cell-inducing capacity of pDCs, which remained significantly higher in comparison to CD1c C DC single cultures (p D 0.0095 for cultures stimulated with CpG-C, p D 0.0004 for cultures stimulated with CpG-C and poly(I:C), p D 0.0043 for cultures stimulated with CpG-C and LPS). Thus, owing to the presence of pDCs, DC co-cultures retain the ability to induce plasma cell differentiation.
Combining CD1c C DCs and pDCs synergistically enhances NK cell cytotoxicity DCs do not only take part in the control of adaptive immunity, but can also interact with innate immune cells such as NK cells. Both IL-12 and IFN-a, cytokines secreted by human blood DCs, are known activators of NK cells. 23 Therefore, we investigated whether a combination of human CD1c C DCs and pDCs is favorable for the induction of NK cell responses. Both CD1c C DCs and pDCs exerted the ability to enhance IFN-g secretion by IL-2-activated NK cells, with poly(I:C)-stimulated CD1c C DCs inducing the highest IFN-g levels (Fig. 5A). Similar trends were seen in NK cell cultures with both DC subsets: conditions in which poly(I:C) was used as DC stimulus induced the highest IFN-g levels, but not to a higher level than adding CD1c C DCs only.
In addition to IFN-g secretion, NK cell cytolytic activity, the hallmark of NK functionality, was determined. Following their activation by DCs, NK-killing of either the NK-sensitive K562 tumor cells or the NK-resistant Daudi tumor cells was determined by flow cytometry. Both CpG-C-stimulated pDCs and poly(I:C)-stimulated CD1c C DCs enhanced the killing of K562 and Daudi tumor cells by autologous NK cells (Fig. 5 B and C). When NK cells were activated by the combination of poly(I:C)stimulated CD1c C DCs and CpG-C-stimulated pDCs, a significant further increase in Daudi cell killing was detected (p D 0.0011 when compared to CpG-C-activated pDCs and p D 0.0115 when compared to poly(I:C)-activated CD1c C DCs; Fig. 5C). To exclude the possibility that this was merely an effect of the combination of CpG-C and poly(I:C), both stimuli were also tested on DC single cultures. When compared to DC single cultures stimulated with both CpG-C and poly(I:C), the increase in Daudi cell killing in the DC co-cultures remained significant (p D 0.0010 when compared to pDCs and p D 0.0002 when compared to CD1c C DCs; Fig. S3).
In an attempt to investigate the mechanism behind the observed synergistic effect on NK-cytotoxicity, other CpG classes were used to stimulate the DCs. Although CpG-A mainly induces the secretion of high IFN-a levels by pDCs, CpG-B mainly elevates the expression of co-stimulatory molecules. 24 As expected, DC co-cultures stimulated with CpG-A secreted more IFN-a than those stimulated with CpG-C, whereas stimulation with CpG-B did not lead to detectable levels of IFN-a (Fig. S4). When CpG-A was used as stimulus in the DC co-cultures instead of CpG-C, higher levels of NK-mediated Daudi cell killing were achieved (p D 0.0248; Fig. 5D). In contrast, the use of CpG-B resulted in less killing, although the difference was not statistically significant (p D 0.1603). Furthermore, addition of recombinant IFN-a to NK cells alone also increased their cytotoxicity (p D 0.0237; Fig. 5E). These results hint toward a role of IFN-a in the synergistic enhancement of NK cell cytotoxicity in co-cultures of human CD1c C DCs and pDCs, which is consistent with the observed increase in IFN-a production in the DC co-cultures. In summary, human blood CD1c C DCs and pDCs could cross activate each other, while combined TLR-stimulation of both subsets induced further changes, including increased CD86 expression on pDCs and enhanced IFN-a production. These crosstalk-mediated effects did not affect the DCs' ability to induce T cell polarization, but did augment NK cell-mediated killing of an NK-resistant tumor cell line. Furthermore, co-application of CD1c C DCs and pDCs allowed for complementation of subset-specific effector functions, such as induction of plasma cell differentiation from B cells and secretion of high levels of IFN-g by PBLs and NK cells.

Discussion
Naturally occurring DCs might become a promising alternative to monocyte-derived DCs for use in DC-based cancer immunotherapy. 3 Our group previously conducted the first clinical trials with ex vivo stimulated and TAA-pulsed primary human DCs in metastatic melanoma patients and demonstrated the safety, feasibility and proof of principle of this approach. 15,16 Overall survival of patients that received a pDC vaccine was greatly increased compared to matched controls and in 7 out of 15 patients an increase in TAA-specific CD8 C T cells was detected after vaccination. 15 Of 14 patients that received a CD1c C DC vaccine, 4 showed long-term progression-free survival, which correlated with multifunctional CD8 C T-cell responses in 3 of these patients. 16 Others have shown that primary DCs might also be used for treatment of cancers other than melanoma. Vaccination with TAA-pulsed CD1c C DCs for advanced stage metastatic prostate cancer was considered feasible and safe. 14 These positive results prompted us to investigate whether a combination of primary DC subsets can further enhance immunological responses.
A drawback of the use of primary DCs for immunotherapy is that the number of DCs that can be isolated from peripheral blood is relatively low. Combining CD1c C DCs and pDCs for vaccination purposes would mean that a larger number of DCs can be adoptively transferred, thereby potentially enhancing immune responses. In addition, CD1c C DCs and pDCs are functionally distinct and may interact differently with adaptive and innate immune cells, suggesting that they might be able to complement each other's effector functions. Furthermore, several in vitro and in vivo studies have demonstrated that direct crosstalk between DC subsets can take place both in the human and murine setting and is mediated by soluble and contactdependent mechanisms. 25,26 Evidence for synergy between mouse DC subsets was found in a study by Lou et al., which showed that immunization of mice with a mixture of mDCs and activated pDCs improved antitumor responses, compared to immunization with either subset alone. 27 In this study, we compared co-cultures of CD1c C DCs and pDCs to single cultures of either DC subset for their phenotype, cytokine secretion and ability to induce effector cell responses. The total amount of DCs was kept the same in DC single cultures and DC co-cultures, to exclude cell number-dependent effects. TLR agonists were added to mature one particular DC subset: TLR-3 agonist poly(I:C) for CD1c C DCs and TLR-9 agonist CpG-C for pDCs. In line with previous findings, 28,29 we found that mature CD1c C DCs can activate immature pDCs and vice versa. This cross activation effect could therefore also play a role during DC vaccination therapy, where infused mature DCs might propagate initial induced immune responses by maturing host DCs.
One of the issues that needs to be considered when using CD1c C DCs and pDCs into a combined vaccine is whether these subsets should be maturated separately or together. Based on the findings in this manuscript, should only one TLR agonist specific to one of the subsets be used, then it is preferable to stimulate the CD1c C DCs and pDCs together, so that the mature subset can mediate cross activation of the other subset. Also in the in vivo study by Lou et al., maturation of pDCs in the presence of mDCs was crucial, as a mixture of mDCs and separately activated pDCs did not result in enhanced antigenspecific CD8 C T cell responses. 27 Since the mDCs in this model were not stimulated by a TLR agonist, the observed synergistic effect can be explained by the cross activation of the immature mDCs by mature pDCs prior to vaccination. In an ideal DC vaccine, however, both subsets should be fully activated, in order to induce maximal immunological responses. Here, we have shown that combined maturation of human blood DC subsets using TLR agonists can also affect their phenotype and functionality. pDCs cultured with CD1c C DCs in the presence of poly(I:C) and CpG-C, upregulated CD86, compared to pDCs receiving similar stimuli in the absence of CD1c C DCs. However, although crosstalk between DC subsets leads to the upregulation of co-stimulatory molecules, it also leads to upregulation of PD-L1. While engagement of co-stimulatory molecules enhances T cell activation, binding of PD-L1 to the PD-1 receptor on T cells provides an inhibitory signal that suppresses T cell activation. 30 Besides changing their expression of surface molecules, co-cultures of CD1c C DCs and pDCs stimulated with poly(I:C) and CpG-C produced more IFN-a than single cultures of stimulated pDCs. The difference was significant even though the absolute number of pDCs, the main producers of IFN-a, in the DC co-cultures was half of the amount present in the DC single cultures. IFN-a can, among other effects, induce Th1 skewing of naive CD4 C T cells, 12 which in turn contributes to the antitumor immune response.
To determine the overall effect of CD1c C DC-pDC crosstalk on T cell responses, we cultured stimulated DCs with naive CD4 C T cells or naive CD8 C T cells. No significant changes in cytokine production by the T cells were observed for the combination of DC subsets compared to the use of single DC subsets. Hence, combining CD1c C DCs with pDCs does not seem to negatively alter their ability to induce CD4 C and CD8 C T cell responses. However, the observed upregulation of PD-L1 in the DC co-cultures, suggests that certain DC vaccines might benefit from the combination with blocking antibodies directed against PD-L1 and/or PD-1. Thus, a combination of DC subsets may synergistically enhance T cell responses, when used in conjunction with these checkpoint inhibitors.
Besides efficiently priming naive T cells, DCs possess the ability to activate other immune effector cells, including B cells and NK cells. NK cells are critical for the antitumor immune response, but the role of B cells in tumor immunity is controversial. [31][32][33][34][35][36][37] Nevertheless, it has recently been reported that antibodies directed against tumor cells can promote antitumor responses by facilitating enhanced uptake of tumor antigens by DCs. 38 In addition, tumor cell-specific antibodies can mediate complement dependent and cell dependent cytotoxicity. 21 Antibodies are mainly produced by plasma cells, which develop from activated B cells. We have shown that pDCs are superior to CD1c C DCs at enhancing plasma cell differentiation from activated B cells. In contrast, CD1c C DCs outperformed pDCs at enhancing IFN-g production in both an MLR assay and in culture with activated autologous NK cells. The advantage of one particular DC subset was retained when the subsets were combined. Thus, human blood mDCs and pDCs can complement each other's effector functions, even when cultured together.
CTLs can detect and lyse transformed cells in an antigendependent manner, by recognizing neo-antigens or overexpressed self-antigens. DC vaccine-based immunotherapy has therefore historically focused at inducing strong CD8 C T cell responses. However, the importance of DC interactions with innate cells for clinical outcome is now becoming more appreciated. 39 NK cells form the body's main line of defense against tumor cells that downregulate MHC class I to evade detection by CTLs. Besides reacting to a loss of inhibitory ligands, NK cells possess activating receptors that recognize stress-induced ligands on tumor cells. NK cells do not necessarily require prior activation to lyse target cells and in our experiments they could efficiently kill K562 tumor cells. Still, addition of poly(I:C)-activated CD1c C DCs, CpG-C-activated pDCs or both to the culture enhanced NK cell cytotoxicity. However, other tumor cells may be more resistant to NK cell-mediated killing. In order to lyse such target cells, NK cell do require prior stimulation. In line with previous findings, 40 we found that NK cells alone only induce minor killing of Daudi tumor cells, but culture of NK cells with poly(I:C)-activated CD1c C DCs or CpG-C-activated pDCs enhanced their killing potential. Moreover, the combination of both DC subsets further increased the killing. This synergistic effect is most likely mediated by the increased amount of IFN-a produced by the DC co-cultures, as replacement of CpG-C by CpG-A to mature DCs further enhanced NK-mediated Daudi cell killing, while use of CpG-B instead lowered the killing. No synergistic increase in NK-mediated killing could be detected in our assays with K562 tumor cells, most likely because optimal killing is already achieved when one of the DC subsets is added.
Our results indicate that combining CD1c C DCs and pDCs for therapeutic purposes can have beneficial effects. Furthermore, stimulation of CD1c C DCs and pDCs together allows for direct crosstalk between the subsets, which can alter their phenotype and functionality. However, using multiple stimuli to facilitate maturation of DC subsets in co-culture may have additional, non-crosstalk-related effects. CD1c C DCs did not optimally mature in the presence of both poly(I:C) and CpG-C: the increase in co-stimulatory molecules was lower than when CpG-C was absent and no detectable levels of IL-12p70 were produced. This can have consequences for their effector function, as DC co-cultures that received poly(I:C) and CpG-C induced lower IFN-g levels in both an MLR assay and in culture with autologous NK cells, than when poly(I:C) and IL-3 were used as stimuli. The inhibitory effect of CpG ODN on poly(I:C)induced DC maturation has been described before for monocyte-derived DCs. 41 We therefore took along LPS as a stimulus for CD1c C when they were subsequently cultured with effector cells. Although CD1c C DCs stimulated with this TLR-4 agonist express high levels of co-stimulatory molecules, they do not secrete detectible levels of IL-12p70. 8 This is likely the reason why LPS-stimulated CD1c C DCs were unable to enhance IFN-g production or tumor cell killing by autologous NK cells. Ideally, the DC subsets should be stimulated in such a way that pDCs secrete high levels of IFN-a and CD1c C DCs secrete high levels of bioactive IL-12. While CpG ODN inhibits poly(I:C)-induced DC maturation, a combination of other TLR agonists may actually enhance DC activation. 42 Alternatively, both DC subsets can be stimulated with a single TLR agonist. Our group recently established a clinical-grade stimulus that can activate both CD1c C DCs and pDCs in the form of stabilized RNA complexes. 43 The use of such a stimulus would be optimal to exploit the advantages of both DC subsets in a clinical setting.
Taken together, a combination of CD1c C DCs and pDCs allows for complementation or even enhancement of each other's effector functions. Hence, combining different human blood DC subsets holds promise for improvement of current DC-based anticancer vaccines.

DC activation
The following TLR ligands were used for DC stimulation: 5 mg/

Cell lines
The human NK-sensitive K562 tumor cell line was obtained from the American Type Culture Collection (ATCC, catalog number: K-562 ATCC CCL-243). The human NK-resistant Daudi tumor cell line was obtained from ATCC (catalog number: Daudi ATCC CCL-213) and kindly provided to us by the laboratory of Prof Kris Thielemans (Free University of Brussels, Brussels, Belgium). Both cell lines were cultured in Iscove's modified Dulbecco's medium (Invitrogen, 12440053) supplemented with 10% fetal bovine serum (Greiner Bio-One) and 1% antibiotic antimycotic (Gibco, 15240062).

Phenotype and cytokine secretion of stimulated DCs
CD1c C DCs and pDCs were incubated overnight at 37 C with different stimuli in a 96-well round-bottom plate, either separate or together at a 1:1 ratio, with each well containing equal cell numbers (50 £ 10 3 cells in 100 mL). After overnight culture, supernatants were taken and cells were stained with the following primary monoclonal antibodies: anti-human leukocyte antigen (

Mixed lymphocyte reaction
CD1c C DCs and pDCs were incubated overnight at 37 C with different stimuli in a 96-well round-bottom plate, either separate or together at a 1:1 ratio, with each well containing equal cell numbers. After overnight culture, DCs were washed with culture medium to remove residual stimuli before subsequent culture with effector cells.
The ability of DCs to induce T cell proliferation and cytokine production was tested in an MLR, using 100 £ 10 3 allogeneic PBLs that were added to the DCs at a 10:1 ratio (PBL:DC) in 200 mL and incubated for 4 d at 37 C. On day 3, supernatants were taken and analyzed for IFN-g (Thermo Fisher Scientific, M700A) by sandwich ELISA. Cellular proliferation was analyzed by measuring overnight [ 3 H]-thymidine (1 mCi/well; MP Biomedicals, 2407005) incorporation on day 4 with a b counter.

Naive CD8 C T cell activation
CD1c C DCs and pDCs were stimulated as above and co-cultured with allogeneic naive CD8 C T cells at a 4:1 ratio (T cell: DC) in the presence of 1 pg/mL SEB. Expression of granzyme B by CD8 C T cells was measured on day 5 by intracellular staining using anti-GrB-PE antibody (Sanquin, M2289). To detect intracellular cytokine production, CD8 C T cells were restimulated on day 6 with PMA/ionomycin in the presence of brefeldin A for 6 h and after fixation and permeabilization, stained with anti-IFN-g-PerCp-Cy5 and anti-TNF-a-PE (Biolegend, 502909). Samples were measured on a CyAn ADP.

Plasma cell differentiation
In contrast to the other effector cell assays, B cells were co-cultured directly with freshly isolated DCs. Autologous B cells were incubated with or without DCs at a 10:1 ratio (B cell:DC) in the presence of different stimuli, for 7 d at 37 C. DCs were CD1c C DCs, pDCs or a combination of both at a 1:1 ratio, with each well containing equal total DC numbers. After 7 d, B cells were stained with CD20-APC (eBioscience, 17-0209-42) and CD38-PE (ImmunoTools, 21270384) or CD38-PerCP-Cy5.5 (BD Biosciences, 551400) and a viability dye (Fixable Viability Dye eFluor 780 or 506; eBioscience, 65-0865-18 or 65-0866-18) and measured on a FACSVerse. Plasma cells were defined as CD38 high CD20 low expressing cells.

NK cell activation and cytotoxicity
CD1c C DCs and pDCs were stimulated as mentioned for the MLR assay, co-cultured with 100 £10 3 autologous NK cells at a 1:1 ratio in 200 mL and incubated for 24 h at 37 C in the presence of 100 U/mL IL-2. NK cells cultured without DCs, but with IL-2 or IL-2 and 1 ng/mL IL-12p70 (BD Biosciences, 554613), were used as controls. Supernatants were taken after 24 h and analyzed for IFN-g by sandwich ELISA.
For cytotoxicity assays, CD1c C DCs and pDCs were stimulated as above, co-cultured with autologous NK cells at a 5:1 ratio (NK:DC) and incubated for 24 h at 37 C. NK cells cultured without DCs and with or without 10 ng/mL IFN-a-2a (Roferon Ò -A; Roche) were used as reference. Target cells (K562 or Daudi) were labeled with PKH67 green fluorescent cell linker (Sigma-Aldrich, MIDI67-1KT) according to manufacturer's instruction and added to the NK cell cultures for the last 4 h of the culture at an NK cell:target cell ratio of 5:1. Target cells cultured without NK cells were used as negative control. Cell death was determined after staining with propidium iodide (PI; Biolegend, 1056-1) and annexin-V-APC (BD Biosciences) in annexin-V binding buffer (BD Biosciences, 550475). Samples were measured on a FACSVerse or FACSAria II (BD Biosciences). Specific killing of target cells was calculated using the formula: %killing D 100 ¡ [(%Annexin¡V ¡ PI ¡ target cells with NK cells/%Annexin¡V ¡ PI ¡ target cells without NK cells) £ 100]. 44 To assess IFN-a secretion by DC co-cultures during the cytotoxicity assay, CD1c C DCs and pDCs were stimulated as above (20 £10 3 cells in 100 mL). After 12 h of overnight culture, supernatants were taken, and DCs were washed with culture medium to remove residual stimuli before subsequent culture of 24 h, after which supernatants were again taken. As a control, supernatants were also taken of DCs cultured for 36 h without washing.

Statistical analysis
GraphPad Prism software (v6.01) was used for statistical analysis. For multiple comparisons, repeated measures one-way ANOVA was performed, followed by post-hoc Dunnett's test when comparisons were made with one group, or Bonferroni's test when comparisons were made with all groups, to detect statistical significant differences between groups. For single comparisons, paired Student's t-tests were performed to detect statistical significant differences between groups. p-Values < 0.05 were considered statistically significant ( Ã p < 0.05; ÃÃ p < 0.01; ÃÃÃ p < 0.001).

Disclosure of potential conflicts of interest
No potential conflicts of interest were disclosed.


extfieldsofstudy:
Biology
Computer Science
Medicine

text:
What time is it? Deep learning approaches for circadian rhythms

Motivation: Circadian rhythms date back to the origins of life, are found in virtually every species and every cell, and play fundamental roles in functions ranging from metabolism to cognition. Modern high-throughput technologies allow the measurement of concentrations of transcripts, metabolites and other species along the circadian cycle creating novel computational challenges and opportunities, including the problems of inferring whether a given species oscillate in circadian fashion or not, and inferring the time at which a set of measurements was taken. Results: We first curate several large synthetic and biological time series datasets containing labels for both periodic and aperiodic signals. We then use deep learning methods to develop and train BIO_CYCLE, a system to robustly estimate which signals are periodic in high-throughput circadian experiments, producing estimates of amplitudes, periods, phases, as well as several statistical significance measures. Using the curated data, BIO_CYCLE is compared to other approaches and shown to achieve state-of-the-art performance across multiple metrics. We then use deep learning methods to develop and train BIO_CLOCK to robustly estimate the time at which a particular single-time-point transcriptomic experiment was carried. In most cases, BIO_CLOCK can reliably predict time, within approximately 1 h, using the expression levels of only a small number of core clock genes. BIO_CLOCK is shown to work reasonably well across tissue types, and often with only small degradation across conditions. BIO_CLOCK is used to annotate most mouse experiments found in the GEO database with an inferred time stamp. Availability and Implementation: All data and software are publicly available on the CircadiOmics web portal: circadiomics.igb.uci.edu/. Contacts: fagostin@uci.edu or pfbaldi@uci.edu Supplementary information: Supplementary data are available at Bioinformatics online.

Introduction
The importance of circadian rhythms cannot be overstated: circadian oscillation have been observed in animals, plants, fungi and cyanobacteria and date back to the very origins of life on Earth. Indeed, some of the most ancient forms of life, such as cyanobacteria, use photosynthesis as their energy source and thus are highly circadian almost by definition. These oscillations play a fundamental role in coordinating the homeostasis and behavior of biological systems, from the metabolic (Eckel-Mahan and Sassone-Corsi, 2009;Froy, 2011;Takahashi et al., 2008;Yoo et al., 2004) to the cognitive levels (Eckel-Mahan et al., 2008;Gerstner et al., 2009). Disruption of circadian rhythms has been directly linked to health problems (Knutsson, 2003;Lamia et al., 2008;Takahashi et al., 2008) ranging from cancer, to insulin resistance, to diabetes, to obesity and to premature ageing (Antunes et al., 2010;Froy, 2010Froy, , 2011Karlsson et al., 2001;Kohsaka et al., 2007;Kondratov et al., 2006;Sharifian et al., 2005;Shi et al., 2013;Turek et al., 2005). At their most fundamental level, these oscillations are molecular in nature, whereby the concentrations of specific molecular species such as transcripts, metabolites and proteins oscillate in the cell with a 24 h periodicity. Modern high-throughput technologies allow largescale measurements of these concentrations along the circadian cycle thus creating new datasets and new computational challenges and opportunities. To mine these new datasets, here we develop and apply machine learning methods to address two questions: (i) which molecular species are periodic? and (ii) what time or phase is associated with high-throughput transcriptomic measurements made at a single timepoint?
At the molecular level, circadian rhythms are in part driven by a genetically encoded, highly conserved, core clock found in nearly every cell based on negative transcription/translation feedback loops, whereby transcription factors drive the expression of their own negative regulators (Partch et al., 2014;Schibler and Sassone-Corsi, 2002), and involving only a dozen genes (Partch et al., 2014;Yan et al., 2008). In the mammalian core clock (Fig. 1), two bHLH transcription factors, CLOCK and BMAL1 heterodimerize and bind to conserved E-box sequences in target gene promoters, thus driving the rhythmic expression of mammalian Period (Per1, Per2 and Per3) and Cryptochrome (Cry1 and Cry2) genes (Stratmann and Schibler, 2006). PER and CRY proteins form a complex that inhibits subsequent CLOCK:BMAL1-mediated gene expression (Brown et al., 2012;Dibner et al., 2010;Partch et al., 2014). The master core clock located in the suprachiasmatic nucleus (SCN) (Moore and Eichler, 1972;Ralph et al., 1990) of the hypothalamus interacts with the peripheral core clocks throughout the body (Takahashi et al., 2008;Yoo et al., 2004).
In contrast to the small size of the core clock, high-throughput transcriptomic (DNA microarrays, RNA-seq) or metabolomic (mass spectrometry) experiments (Andrews et al., 2010;Eckel-Mahan et al., 2012Hughes et al., 2009;Masri et al., 2014b;Miller et al., 2007;Panda et al., 2002;P.Tognini et al., in preparation), have revealed that a much larger fraction, typically on the order of 10%, of all transcripts or metabolites in the cell are oscillating in a circadian manner. Furthermore, the oscillating transcripts and metabolites differ by cell, tissue type, or condition (Panda et al., 2002;Storch et al., 2002;Yan et al., 2008). Genetic, epigenetic and environmental perturbations-such as a change in diet-can lead to cellular reprogramming and profoundly influence which species are oscillating in a given cell or tissue (Bellet et al., 2013;Dyar et al., 2014;Eckel-Mahan et al., 2012Masri et al., 2013Masri et al., , 2014a. When results are aggregated across tissues and conditions, a very large fraction, often exceeding 50% and possibly approaching 100%, of all transcripts is capable of circadian oscillations under at least one set of conditions, as shown in plants (Covington et al., 2008;Harmer et al., 2000), cyanobacteria and algae (Monnier et al., 2010;Vijayan et al., 2009) and mouse (Patel et al., 2015;Zhang et al., 2014).
In a typical circadian experiment, high-throughput omic measurements are taken at multiple timepoints along the circadian cycle under both control and treated conditions. Thus the first fundamental problem that arises in the analysis of such data is the problem of detecting periodicity, in particular circadian periodicity, in these time series. The problem of detecting periodic patterns in time series is of course not new. However, in the cases considered here the problem is particularly challenging for several reasons, including: (i) the sparsity of the measurements (the experiments are costly and thus data may be collected for instance only every 4 h); (ii) the noise in the measurements and the well known biological variability; (iii) the related issue of small sample sizes (e.g. n ¼ 3); (iv) the issue of missing data; (v) the issue of uneven sampling in time; and (vi) the large number of measurements (e.g. 20 000 transcripts) and the associated multiple-hypothesis testing problem. Here we develop and apply deep learning methods for robustly assessing periodicity in high-throughput circadian experiments, and systematically compare the deep learning approach to the previous, non-machine learning, approaches (Glynn et al., 2006;Hughes et al., 2010;Yang and Su, 2010). While this is useful for circadian experiments, the vast majority of all high-throughput expression experiments have been carried, and continue to be carried, at single timepoints. This can be problematic for many applications, including applications to precision medicine, precisely because circadian variations are ignored creating possible confounding factors. This raises the second problem of developing methods that can robustly infer the approximate time at which a single-time high-throughput expression measurement was taken. Such methods could be used to retrospectively infer a time stamp for any expression dataset, in particular to improve the annotations of all the datasets contained in large gene expression repositories, such as the Gene Expression Omnibus (GEO) (Edgar et al., 2002), and improve the quality of all the downstream inferences that can be made from this wealth of data. There may be other applications of such a method, for instance in forensic sciences, to help infer a time of death. In any case, to address the second problem we also develop and apply deep learning methods to robustly infer time or phase for single-time high-throughput gene expression measurements.

Periodicity inference from time series measurements
To train and evaluate the deep learning methods, we curate BioCycle, the largest dataset including both synthetic and real-world biological time series, and both periodic and aperiodic signals. While the main goal here is to create methods to analyze real-world biological data, relying only on biological data to determine the effectiveness of a method is not sufficient because there are not many biological samples which have been definitively labeled as being periodic or aperiodic. Even when one can be confident that a signal is periodic, it can be difficult to determine the true period, phase and amplitude of that signal. Therefore, we rely also on synthetic data to provide us with signals that we can say are definitely periodic or aperiodic, and whose attributes-such as period, amplitude, and phase-can be controlled and are known. Furthermore, previous approaches were developed using synthetic data and thus the same synthetic data must be used to make fair comparisons.

Synthetic data
We first curate a comprehensive synthetic dataset BioCycle Synth, which includes all previously defined synthetic signals found in JTK_Cycle (Hughes et al., 2010) and ARSER (Yang and Su, 2010), but also contains new signals. BioCycle Synth is in turn a collection of two different types of datasets: a dataset in which signals are constructed using mathematical formulas (BioCycle Form) , and a dataset in which signals are generated from a Gaussian process (Rasmussen, 2004) (BioCycle Gauss) . In previous work, synthetic data was generated with carefully constructed formulas to try to mimic periodic signals found in real-world data (see below). While this gives one a lot of control over the data, it can create signals that are too contrived and therefore not representative of real-world biological variations. In addition, the noise added at each timepoint is independent of the other timepoints, which may not be the case in real-world data. The BioCycle Gauss dataset uses Gaussian processes to generate the data and address these problems.
The datasets used in JTK_Cycle contain the following types of formulas or signals: cosine, cosine with outlier timepoints and white noise. The ARSER dataset contains cosine, damped cosine with an exponential trend, white noise and an auto-regressive process of order 1 (AR(1)). In addition to all the aforementioned signals, BioCycle Form contains also 9 additional kinds of signals: combined cosines (cosine2), cosine peaked, square wave, triangle wave, cosine with a linear trend, cosine with an exponential trend, cosine multiplied by an exponential, flat and linear signals (many of which can be found in Deckard et al., 2013). Figure 2 shows an example of each type of signal found in the BioCycle Form dataset. For clarity, the periodic signals are shown without noise. Signals in the BioCycle Form dataset have an additional random offset chosen uniformly between À200 and 200, random amplitudes chosen uniformly between 1 and 100, signal to noise ratios (SNRs) of 1-5, random phases chosen uniformly between 0 and 2p, and periods between 20 and 28. [Data for the second (12 h) and third (8 h) harmonics, which are found in biological data, are also generated (Supplementary Information)]. At each timepoint sample, zero mean Gaussian noise is added with the proper SNR variance.
The BioCycle Gauss dataset is obtained from a Gaussian process. The value of the covariance matrix corresponding to the timepoints x and x 0 is determined by a kernel function kðx; x 0 Þ. Equation 1 is the kernel function used to generate the periodic signals, and Equation 2 is the kernel function used to generate the aperiodic signals in BioCycle Gauss .
The parameter l controls how strong the covariance is between two different timepoints, r controls how noisy the synthetic data is, and b can add a non-stationary, linear, trend to the signals (Duvenaud, 2014). The parameter p in equation 1 is the period of the signal. To generate the data in BioCycle Gauss, the values of l, r, b, p, as well as the offset and the scale are varied, in a way similar to the data in BioCycle Form . Examples of signals from the BioCycle Gauss dataset are given in Figure 3.
JTK_Cycle analyzes synthetic signals sampled over 48 h with a sampling frequency of 1 and 4 h. ARSER analyzes synthetic signals sampled over 44 h with a sampling frequency of 4 h. BioCycle analyzes synthetic signals sampled over 24 and 48 h. Signals sampled over 24 h have a sampling frequency of 4, 6 and an uneven sampling at timepoints 0, 5, 9, 14, 19 and 24. Signals sampled over 48 h have sampling frequencies of 4, 8 and an uneven sampling at timepoints 0,4,8,13,20,24,30,36,43. The sampling frequencies in these datasets are intentionally sparse to mimic the sparse temporal sampling of real-world high-throughput data. The number of synthetic signals at each sampling frequency is 1024 for JTK_Cycle, 20 000 for ARSER and 40 000 for BioCycle Synth . Finally, each signal in BioCycle Synth has three replicates, obtained by adding random Gaussian noise to the signal, to mimic typical biological experiments.

Biological data
The performance of any circadian rhythm detection method requires extensive validation on biological datasets. In previous work, due to the aforementioned difficulty of not having ground truth labels, the biological signals detected as being periodic had to be inspected by hand, or loosely assessed by comparison to other methods (Straume, 2004). In addition to the scaling problems associated with manual inspection, this approach did not allow the computation of precise classification metrics (Baldi et al., 2000), such as the AUC-the Area Under the Receive Operating Characteristic (ROC) Curve. The repository of circadian data hosted on CircadiOmics (Patel et al., 2012) includes over 30 high-throughput circadian transcriptomic studies, as well as several circadian high-throughput metabolomic  studies, that provide extensive coverage of different tissues and experimental conditions. From the CircadiOmics data, a high-quality biological dataset BioCycle Real is created with periodic/aperiodic labels.
To curate BioCycle Real, we start from 36 circadian microarray or RNA-seq transcriptome datasets, 32 of which are currently publicly available from the CircadiOmics web portal (28 of these are also available from CircaDB (Pizarro et al., 2012)). Five datasets are from ongoing studies and will be added to CircadiOmics upon completion. All included datasets correspond to experiments carried out in mice, with the exception of one dataset corresponding to measurements taken in Arabidopsis Thaliana. BioCycle Real comprises experiments carried over a: 24-h period with a 4 h sampling rate; 48-h period with a 2 h sampling rate; and 48-h period with a 1 h sampling rate.
To extract from this set a high-quality subset of periodic time series, we focus on the time series associated with the core clock genes (Fig. 1) in the control experiments. These gene include Clock, Per1, Per2, Per3, Cyr1, Cry2, Nr1d1, Nr1d2, Bhlhe40, Bhlhe41, Dbp, Npas2 and Tel (Harmer et al., 2001) for mouse, and the corresponding orthologs in Arabidopsis (Harmer et al., 2000). Arabdiposis orthologs were obtained from Affymetrix NetAffx probesets and annotations (Liu et al., 2003). These core gene time series were further inspected manually to finally yield a set of 739 high-quality periodic signals. To extract a high-quality biological aperiodic dataset, we start from the same body of data. To identify transcripts unlikely to be periodic, we select the transcripts classified as aperiodic consistently by all three programs JTK_Cycle, ARSER and Lomb-Scargle with an associated P-value of 0.95. After further manual inspection, this yields a set of 18 094 aperiodic signals. Examples of signals taken at random from the BioCycle Real are shown in Figure 4.

Time inference from single timepoint measurements
To estimate the time associated with a transcriptomic experiment conducted at a single timepoint, we curate the BioClock dataset starting from the same data in CircadiOmics, focusing on mouse data only for which we have enough training data. While in principle inference of the time can be done using the level of expression of all the genes, exploratory feature selection and data reduction experiments (not shown) show that in most cases it is sufficient to focus on the set of core clock genes, or even a subset (see Section 4). Thus the reduced BioClock dataset contains microarray and RNA-Seq single time measurements for each gene transcript in the core clock with the associated timepoint. The BioClock dataset is organized by tissue and condition. Tissues include liver, kidney, heart, colon, glands (pituitary, adrenal), skeletal muscle, bone, white fat and brown fat. Brain specific tissues include SCN (Suprachiasmatic nucleus), hippocampus, hypothalamus and cerebellum. There are also several cell-specific datasets including mouse fibroblasts and macrophages. All the datasets in BioClock contain both control and treatment conditions. There is great variability among the treatment conditions (e.g. Eckel-Mahan et al., 2013;Masri et al., 2014a), varying from gene knock out and knock down (SIRT1 and SIRT6), to changes in diet (high fat, ketogenic), to diseases (epilepsy). It is important to be able to assess the ability of a system to predict time across tissues and conditions.

Methods
We experimented with several machine learning approaches for the two main problems considered here. In general, the best results were obtained using neural networks. This is perhaps not too surprising since it is well known that neural networks have universal approximation properties and deep learning has led to state-of-the art performance, not only in several areas of engineering (e.g. computer vision, speech recognition, natural language processing, robotics) (Hannun et al., 2014;Lenz et al., 2015;Szegedy et al., 2014), but also in the natural sciences Di Lena et al., 2012;Lusci et al., 2013;Quang et al., 2015). Thus here we focus exclusively on deep learning approaches to build two systems, BIO_CYCLE and BIO_CLOCK, to address the two main problems. However, we add comparisons to k-nearest neighbors and Gaussian processes in the Supplementary Information.

Periodicity inference from time series measurements
3.1.1 Classifying between periodic and aperiodic signals To classify signals as periodic or aperiodic, we train deep neural networks (DNNs) using standard gradient descent with momentum (Rumelhart et al., 1988;Sutskever et al., 2013). We train separate networks for data sampled over 24 and 48 h. The input to these networks are the expression time-series levels of the corresponding gene (or metabolite). The output is computed by a single logistic unit trained to be 1 when the signal is periodic and 0 otherwise, with relative entropy error function. We experimented with many hyperparameters and learning schedules. In the results reported, the learning rate starts at 0.01, and decays exponentially according to 0:1 1:000087 t , where t is the iteration number. The training set consists of 1 million examples, a size sufficient to avoid overfitting. The DNN uses a mini-batch size of 100 and is trained for 50 000 iterations. Use of dropout (Baldi and Sadowski, 2014;Srivastava et al., 2014), or other forms of regularization, leads to no tangible improvements. The best performing DNN found (Fig. 5(a)) has 3 hidden layers of size 100. We are able to obtain very good results by training BIO_CYCLE on synthetic data alone and report test results obtained on BioCycle Form, BioCycle Gauss and BioCycle Real . 3.1.2 Estimating the period In a way similar to how we train DNNs to classify between periodic and aperiodic signals, we can also train DNNs to estimate the period of a signal classified as periodic. During training, only periodic time series are used as input to train these regression DNNs. The output of the DNNs are implemented using a linear unit and produce an estimated value for the period. The error function is the squared error between the output of the network and the true period of the signal, which is known in advance with synthetic data. Except for the difference in the output unit, we use the same DNNs architectures and hyperparameters as for the previous classification problem.

Estimating the phase and the lag
After the period p, we estimate the phase / of a signal s by finding the value / that maximizes the following expression: P t2T cosð 2pt p þ /Þs½t, where T is the set of all timepoints. Given /, the lag (i.e. at what time the periodic pattern starts) is given by /p 2p .

Estimating the amplitude
After the phase /, we estimate the amplitude a by first removing any linear trend and then comparing the variance of the signal to the variance of a cosine signal with parameters /, p and amplitude 1. The formula is shown in Equation 3, where l s ¼ 1 jTj P t2T s½t and l c ¼ 1 We cannot claim this approach is new, however, we have not seen it in previous literature. An alternative is to measure the amplitude on the smoothed time series.

Calculating P-values and q-values
To calculate P-values, the distribution of the null hypothesis must first be obtained. To do this, N aperiodic signals are generated from one of the two BioCycle Synth datasets. Then we calculate the N output values V(i) (i ¼ 1; . . . ; N) of the DNN on these aperiodic signals. The P-value for a new signal s with output value V is now 1 N X N i¼1 1ðV > VðiÞÞ, where 1 is the indicator function. This equation provides an empirical frequency estimate for the probability of obtaining an output of size V or greater, assuming that the signal s comes from the null distribution (the distribution of aperiodic signals). Therefore, the smaller the P-value, the more likely it is that s is periodic. The q-values are obtained through the Benjamini and Hochberg procedure (Benjamini and Hochberg, 1995). We also compute a posterior probability of periodic expression (PPPE) using the method described in Allison et al. (2002), which models the distribution of P-values as a mixture of beta distributions.

Time inference from single timepoint measurements
For this task, different machine learning methods were investigated, including simple linear regression, k-nearest neighbors, decision trees, shallow learning and deep learning, including unsupervised compressive autoencoders (Baldi, 2012) with two coupled phase (cosine/sine) units in the bottleneck layer (below and Supplementary  Information). Supervised deep learning methods give the best results and are used in the final BIO_CLOCK system. The output of the DNNs is implemented using two coupled output units, representing the cosine and the sine of the phase angle ( Fig. 5(b)). If the total weighted inputs into these two units are S 1 and S 2 respectively, then the values of the two outputs units are given by: . These are then automatically converted into a time (ZT). In order to better assess the effect of having data from different tissues, we experiment with both training specialized predictors trained on data originating from a single tissue, as well as predictors trained on data from all tissues. The final general-purpose predictor corresponds to a DNN trained on all the data. In each one of these experiments, we use 5-fold cross validation on the corresponding subset of the BioClock dataset, using architectures with 2 to 9 layers, and 100 to 600 units, to select the best network. A learning rate of 0.1 is typically used, with an exponential decay according to 0:1 1:002 t . A visualization of the DNN is provided in Figure 5(b).

Data normalization
For both the periodicity and time inference problems, training and testing examples are normalized to have a mean of zero and a standard deviation of one.

Software and run time
Downloadable software is currently written in R and Python and is intended to be easy for biologists to use. While exploring different models both Pylearn2 (Goodfellow et al., 2013) and Caffe (Jia et al., 2014) were used. The DNNs typically take hours for training but, once trained, can process a real-world dataset ($20 000 time series) in about one minute, both run times corresponding to a single CPU.

Results
In all the tables, the best results are shown in bold. To determine if the BIO_CYCLE results are significantly different from other methods, the testing set is randomly split into 10 BIO_CYCLE. The output is either the binary periodic/aperiodic classification, or the regression estimate of the period of the signal.

Periodic/aperiodic classification
BIO_CLOCK. The outputs are the cosine and sine of the phase angle associated with the expression measurement of the core clock genes. equal-size, non-overlapping, subsets and the results from each subset are obtained. Then, a Student's t-test is performed between the results of the best of the two DNNs and the best of the previously existing methods. Finally, the P-value from that test is obtained to assess if the result differences are statistically significant. Small Pvalues (such as 0.05 and below) indicate that there is a significant difference between the methods. The P-values from the t-tests are shown in the rightmost column in all the tables. The results focus on periodic signals with periods around 24 h, the most common case, however periods of 12 and 8 h, corresponding to the second and third harmonics, are analyzed in the Supplementary Information. In the tables, the datasets BioCycle Form, BioCycle Gauss and BioCycle Real are referred to as BC F, BC G and BC R, respectively.

Synthetic data (BioCycle Synth)
Results for the area under the receiver operating characteristic curve (AUC) for the task of classifying signals as periodic or aperiodic are shown in Table 1, and the ROC curves computed on BioCycle Form are shown in Figure 6. The DNN F label corresponds to the DNN that has been trained on the BioCycle Form data, and the DNN G label corresponds to the DNN that has been trained on the BioCycle Gauss data. The ROC curves computed on BioCycle Gauss are similar (not shown). The results from Table 1 show that the DNN method has better AUC than all the other published methods on the BioCycle Form and BioCycle Gauss datasets. Though the DNN does better when tested on data from the same distribution as it was trained on, it still outperforms all the other previous methods, regardless of which data it is trained on. A plot showing how the signal to noise ratio (SNR) affects performance is shown in Figure 7. This plot cannot be done for the BioCycle Gauss dataset, since in this case the exact SNR is not known. The DNN outperforms all the other published methods at all SNRs.

Biological data (BioCycle Real )
The performance on the biological dataset is shown in Table 2. Although the ARSER, LS and JTK_Cycle methods achieve good performance on the aperiodic data, as can be expected since they were used to label the aperiodic data, the DNN method remains very competitive, often outperforming at least one of the other published methods.

Evaluation of P-value cutoffs
To investigate if the p-values obtained by BIO_CYCLE are reasonable, the accuracy of the periodic/aperiodic classification at different p-value cutoffs is evaluated. In addition to a p-value, BIO_CYCLE produces a binary classification. If the output of the DNN is greater than 0.5 the signals is labeled as periodic, otherwise, it is labeled as aperiodic. The accuracy using this binary classification is also evaluated. Results are shown in Figure 8. The vertical dashed line corresponds to a common p-value cutoff of 0.05. However, a proper p-value does not guarantee that the best accuracy will be at the cutoff of 0.05. Results show that BIO_CYCLE has the highest potential accuracy. It also has the best accuracy at 0.05 for 2 out of the 4 plots. In addition, the binary classification of BIO_CYCLE is almost always better than the accuracy of all the other methods at any p-value cutoff. Histograms showing the p-values can be found in the Supplementary Information.

Period, lag and amplitude estimation
The metric to determine how well each method estimates the period, lag and amplitude is given by the coefficient of determination R 2 .  The line y ¼ x corresponds to perfect prediction. In this case, y is the estimated value given by the method and x is the true value. R 2 measures how well the line y ¼ x fits the points that correspond to the true value versus the estimated value. Perfect prediction corresponds to y ¼ x and corresponds to R 2 ¼ 1. The results for estimating the period, lag and amplitude are shown in Tables 3-5, respectively. For the BioCycle Gauss dataset we cannot control or know the exact lag or amplitude, so there are no results for BioCycle Gauss in Tables 4 and 5. These tables tell a similar story as Table 1. The DNN outperforms the other methods in most of the categories. Even when the DNN is tested on data associated with a distribution that is different from the distribution of its training set, in the majority of the cases it gives superior performance compared to ARSER, LS and JTK_Cycle.

Missing replicates and missing data
In gene expression experiments, replicate measurements can be missing. To investigate how missing replicates affect performance, the BioCycle Form dataset which has three replicates for each timepoint was used to assess performance with zero replicates removed at each timepoint, one replicate removed at each timepoint and two replicates removed at each timepoint. The results are shown in Figure 9 and show that JTK_Cycle is significantly affected in a negative way by missing replicates, while the performance of all the other methods degrades gracefully with the number of missing replicates, and minimally compared to JTK_Cycle. Missing data (timepoints at which there are no replicates) is also handled gracefully by BIO_CYCLE, while it is not handled at all by some of the other methods (not shown).

Time inference from single timepoint measurements
4.4.1 Overall performance BIO_CLOCK is trained using 16 core clock genes: Arntl, Per1, Per2, Per3, Cyr1, Cry2, Nr1d1, Nr1d2, Bhlhe40, Bhlhe41, Dbp, Npas2, Tef, Fmo2, Lonrf3 and Tsc22d3. When trained and tested on all the data, using 70% of the data for training and the remaining 30% for testing, it accurately predicts the time of the experiment with a mean absolute error of 1.22 h (less than 75 min) ( Table 6). We experimented also with training BIO_CLOCK with an even smaller   Table 6 shows the mean absolute errors obtained when training BIO_CLOCK on data from certain organs/tissues and testing it on data from a different set of organs/tissues. All the data is from mice and under WT condition. The only datasets for which we have enough data for training correspond to liver and brain (when aggregating all the corresponding datasets). We form two additional sets (Set 1 and Set 2) by combining data from other organs. The first corresponds to combined data from the adrenal gland, fat, gut, kidney, lung and muscle (Set 1). The second corresponds to combined data from the aorta, colon, fibroblast, heart, macrophages and pituitary gland (Set 2). Finally, all of the aforementioned data is combined to form a bigger dataset (All). In all the experiments reported in Table 6, the data are split using a 70/30 training/test ratio, and tests sets never overlap with any of the corresponding training sets. The DNNs perform best when trained and tested on the same organ/tissue or sets of organ/tissues or when trained on all the organs/tissues. The DNNs perform significantly worse when trained and tested on data with diverging origins. However, in all cases, the DNN trained on the combined dataset does almost as well as, or better than, the corresponding specialized DNN.

Training and testing on different conditions
The collected data also includes data from mice under experimental conditions. The experimental conditions include high-fat and ketogenic diets, epilepsy and SIRT1 and SIRT6 knockouts. This dataset is too small to build a training and testing set. However, one can test the BIO_CLOCK DNN trained on the combined mice organs under normal conditions on this dataset. This experiment yields a mean absolute error of 2.57 h.

Annotation of the GEO database
Finally, we extract all the mouse gene expression experiments contained in the GEO repository (Edgar et al., 2002) and run BIO_CLOCK on them. A file containing all the corresponding imputed times is available from the CircadiOmics web portal.

Conclusion
Deep learning methods can be applied to high-throughput circadian data to address important challenges in circadian biology. In particular, we have developed BIO_CYCLE to detect molecular species that oscillate in high-throughput circadian experiments and extract the characteristics of these oscillations. Remarkably, BIO_CYCLE can be trained with large quantities of synthetic data preventing any kind of overfitting. We have also developed BIO_CLOCK to infer the time at which a transcriptomic sample was collected from the level of expression of a small number of core clock genes. Both methods will be improved as more data becomes available and, more generally, deep learning methods are likely to be useful to address several other related circadian problems, such as analyzing periodicity in high-throughput circadian proteomic data, or inferring sample time in different species. In particular, developing methods for annotating the time of all the human gene expression experiments, contained in GEO, and other similar repositories, would be valuable. Such annotations could be important for improving the interpretation of both old and new data and discovering circadian driven effects that may be important in precision medicine and other applications, for instance to help determine the optimal time for administering certain drugs.