Singularity response reveals entrainment properties in mammalian circadian clock

Entrainment is characterized by phase response curves (PRCs), which provide a summary of responses to perturbations at each circadian phase. The synchronization of mammalian circadian clocks is accomplished through the receipt of a variety of inputs from both internal and external time cues. A comprehensive comparison of PRCs for various stimuli in each tissue is required. Herein, we demonstrate that PRCs in mammalian cells can be characterized using a recently developed estimation method based on singularity response (SR), which represents the response of desynchronized cellular clocks. We confirmed that PRCs can be reconstructed using single SR measurements and quantified response properties for various stimuli in several cell lines. SR analysis reveals that the phase and amplitude after resetting are distinguishable among stimuli. SRs in tissue slice cultures reveal tissue-specific entrainment properties. These results demonstrate that SRs can be employed to unveil entrainment mechanisms with diverse stimuli in multiscale mammalian clocks.

-The results for MEF cells do not include comparison with the classical PCR, but rather focus on the dependence of the phase response on the stimulus concentration, in particular for Dexamethasone (Fig. 5e,f). This is an interesting application of the SR method, where the SR phase and amplitude parameters are used to characterize phase response as a function of DEX concentration. Since the comparison is among SR parameters only, some meaningful information can be obtained (see also my next comment).
-Regarding the SR method itself, it can be shown that the SR phase parameter is equivalent to the phase parameter of the PCR method, i.e. theta(PCR)=theta(SR). This is, however, not true for the amplitude parameter, since a relation of the form R(PCR)=beta*R(SR) needs to be computed for each cell type and the factor beta is different from cell line to cell line, and probably different for each experiment. Therefore, the interpretation of the amplitude parameter for SR can be a problem when recovering PCR and would require further study.
Nevertheless, in the case of comparison among SR experiments only (cf. the previous comment), the values R(SR) retain their meaning and can be used to compare and characterize results from similar experiments, such as the dependence of phase response on the concentration of the stimulus (see Fig.  5e,f). Indeed, as the authors suggest in the Discussion, the SR method could be more suitable for a preliminary, and faster, phase response analysis in experiments such as high throughput screening, where large amounts of tests need to performed. The SR would enable a selection of the best cases, which would then be followed by a more detailed classical PCR. The authors apply their recently published (Masuda et al. Nat. Commun. 2021) singularity response (SR) method for inferring the phase response curve (PRC) from a single experiment to mammalian rather than plant circadian clocks. The method is successfully applied to experimental data in cultured mouse and rat fibroblast cells, as well as several mouse tissue slices for many different stimuli including temperature, PMA, DEX, Forskolin, sodium bicarbonate, LiCl, and CoCl2. This variety of cell types and stimuli provide compelling evidence that the method is broadly useful. This is exciting since it promises to reveal the very insightful PRC with a single experiment from a desynchronized state rather than numerous experiments from a synchronized state for many different phases of the input.
After reading the methods section multiple times and consulting the previous paper on the SR method (Masuda et al. Nat. Commun. 2021), the SR method remains opaque. As I understand it, the method steps are 1) measurements of the individual cell rhythms give estimates of the order parameter R and phase Omega; 2) scale R by beta to account for 'experimental condition' (Omega is left unchanged); 3) assuming the population is desynchronized, invert Eq.(4) to find the PRC g(theta) or more specifically the parameters (a, alpha) in the underlying phase oscillator model that gives rise to the PRC. So that this method can be useful, it is imperative that the exposition of the method be improved greatly. It should be possible for a reasonably competent reader to implement the method for their own data. Here are some specific reasons why I find the method difficult to understand: a) The parameters omega and delta t in Eq.(4) are not defined until much later (after Eq.(12)) in a different context. b) Possibilities for the PRC g(theta) are given in Eq.(5) and Eq.(6), but those are for fitting to PRC data. The g(theta) resulting from the SR method is implicit in the phase oscillator model Eqs.(9) -(12). Perhaps give the expression for the PRC referenced in (17). c) All of the validation results of the method come before the method is adequately explained in the methods section. Perhaps the introduction could include a better summary of the details of the method or an example using synthetic data (generate 'data' from the phase oscillator model and then use the method to infer the known PRC). d) Use the index j instead of i to avoid confusion with the imaginary unit i.
e) The SR method is better explained in Masuda et al. Nat. Commun. 2021. For completeness, the details there should be duplicated here.
f) The details of the inversion process (step 3 above) are not described. This is a non-trivial root finding problem especially since the integral in Eq.(4) is a smoothing operation. I am surprised that the authors do not report experiencing any numerical difficulties (and how they were overcome) with this inversion.
In short, my primary suggestion is that the clarity of the description of the SR method be improved. The paper would also benefit from more investigation of the mechanism behind and the limits of the SR method. Why does it work? When does it not work?
Minor suggestions: 1) Change "mouse" to "mammalian" in the title since the approach was applied to multiple mammalian clocks.

Responses to reviewer's comments.
We thank all the reviewers for their constructive and valuable comments, which have helped us improve our manuscript significantly. We have carefully addressed each comment and revised the manuscript accordingly. The following are the authors' pointby-point responses to the reviewers' comments. We have earnestly tried to address all the concerns raised by the reviewers.

Reviewer #1 (Remarks to the Author):
This manuscript illustrates a new method for measuring circadian clock phase response curves (PRC), which was previously demonstrated in plants (the first author of the current manuscript is also the first author in the previous paper, [17] of the references). The method is based on the Singularity Response (SR) and it consists in applying a stimulus to a de-synchronized cell population. From the phase and amplitude of the population after stimulation, a mathematical model is then applied to recover the PRC parameters.
In principle, the SR method is less costly and of more practical application than the classical PRC method, but it doesn't recover as many details.
The current work tests and illustrates this method for mammalian (mouse) cells, using several cell lines (NIH3T3, MEF, Rat1), as well as tissue slices; several different stimulus according to the cell line; and different clock reporters, also according to cell lines (Reverb, Per, Bmal1). The results are very interesting and comprehensive taking into account the range of cells, stimulus, and reporters tested. In a general way, the new method is in a fair-to-good correspondence with the more classical PRC methods, but when it comes to closer comparison of detailed properties of the SR and PCR, further work and studies are needed. Throughout the manuscript, some points could be further developed: -In the validation for NIH3T3 cells, the authors perform tests with "small populations of 5 cells". It is not clear to me what is the goal of this test. It seems to me that 5 elements is too small for a population, I would expect groups of 20-30 cells in order to obtain reasonable averages. And how do you choose these "cell populations" (are they near neighbors?). On the other hand, how do "populations of 5 cells" compare with the population sizes in the other experiments, in other cell lines? Is it reasonable to assume that these small population results can be extrapolated to "real" populations?
Thank you for your insightful comment. In the original study, we created a small population containing five cells because we wanted to compare the phase response between highly synchronized populations and desynchronized populations (original Figure 1e). When we create a large population, the synchronization rate should be close to the average synchronization rate of all cells, which is extremely low in this case. Thus, we intended to create relatively small populations. However, for the revised manuscript, we analyzed again the phase response in populations composed of 3, 5, 10, and 20 cells to confirm our hypothesis (New Supplementary Figure 4a). We concluded that the amplitude dependency of the phase responses at the cellular and population levels was also close to the mathematical model, as depicted in Fig. 1b. We randomly selected cells and did not consider the spatial relationship among these cells when we created a small population (Fig. 1e, S4b). However, it is also possible to intentionally select cells in close phases to construct a highly synchronized population.
In this case, PRCs are also close to the theoretical model (New Supplementary Figure 4b,Page 6,. A detailed method has been added to the relevant section of the revised manuscript. -The results for Rat-1 cells are not always clear to understand and/or explain. First, the phases computed through SR seem to always be underestimated (all dots above the y=x in Fig 4d). Is this a specific problem for these type of cells?
In the original study, the phase parameter of SR for melatonin was not consistent with that of PRC, while other stimulation showed a similar phase relationship. We repeated the comparison of SR and PRC parameters for melatonin with strictly managed experimental conditions (almost same passage numbers, same lot of culture medium etc.). We found that the response of melatonin was weak, leading to inconsistency in the phase of SR and PRC. As we described, the SR method is useful when the stimulation is strong.
We also measured the PRC parameter for PMA again to confirm whether the underestimation of the SR phase is specific to Rat-1 cells. However, we obtained a similar result, in which the SR phase was slightly delayed compared to the PRC phase (Fig. 4d).
It is possible that the SR phase tends to be underestimated, although the SR phase is almost compatible with the PRC phase.
Second, the PRC and SR parameters for Forskolin and PMA are very close to each other --as indicated by the corresponding dots in Figs. 4c,d.
However, the reconstructed PRC curve for PMA fit the points very badly, in contrast to the Forskolin curve (Fig. 4e). Moreover, the PMA and Forskolin stimulation curves are similar in amplitudes (both high ater stimulation), which seems to indicate that PMA should also provide good results, as was indeed the case in NIH3T3 cells. Can these divergences be explained in some way?
We again measured the PRC for PMA with more strictly managed experimental conditions (almost same passage numbers and the same lot of culture medium as SR measurement, etc.). We found that predicted PRC using the SR method was similar to measured PRC. We speculated that the inconsistency in the two PRCs in the original study was caused by the small number of samples for measured PRC or differences in experimental conditions such as the passage numbers.
-The results for MEF cells do not include comparison with the classical PCR, but rather focus on the dependence of the phase response on the stimulus concentration, in particular for Dexamethasone (Fig. 5e,f). This is an interesting application of the SR method, where the SR phase and amplitude parameters are used to characterize phase response as a function of DEX concentration. Since the comparison is among SR parameters only, some meaningful information can be obtained (see also my next comment).

As per your suggestion, we have added the results of PRC measurements in MEFs treated
with DEX (New Supplementary Fig. 12), confirming the reliability of the SR method, to the revised manuscript.
-Regarding the SR method itself, it can be shown that the SR phase parameter is equivalent to the phase parameter of the PCR method, i.e. theta(PCR)=theta(SR). This is, however, not true for the amplitude parameter, since a relation of the form R(PCR)=beta*R(SR) needs to be computed for each cell type and the factor beta is different from cell line to cell line, and probably different for each experiment. Therefore, the interpretation of the amplitude parameter for SR can be a problem when recovering PCR and would require further study.
Nevertheless, in the case of comparison among SR experiments only (cf. the previous comment), the values R(SR) retain their meaning and can be used to compare and characterize results from similar experiments, such as the dependence of phase response on the concentration of the stimulus (see Fig. 5e,f). Indeed, as the authors suggest in the Discussion, the SR method could be more suitable for a preliminary, and faster, phase response analysis in experiments such as high throughput screening, where large amounts of tests need to performed. The SR would enable a selection of the best cases, which would then be followed by a more detailed classical PCR.
As we described in the supplementary text, beta can be estimated using the SR parameters when the stimuli are strong enough to produce type 0 resetting (i.e. R(PRC) is 1). Herein, we confirmed that 100 nM DEX treatment actually resulted in the type 0 resetting using the conventional method (New Supplementary Fig. 12) and that beta estimated using only the SR parameter was accurate. We apologize for the confusion. We have edited the captions according to your comments. Figs 2g and 4c: give a better idea of what the "approximate line" is, or at least refer to eq.
We have edited the description.

Line 562 (Methods): Should be Fig 2g? (and not 2f)
We have corrected the manuscript as indicated.

Reviewer #2 (Remarks to the Author):
The authors apply their recently published (Masuda et al. Nat. Commun. 2021) singularity response (SR) method for inferring the phase response curve (PRC) from a single experiment to mammalian rather than plant circadian clocks. The method is successfully applied to experimental data in cultured mouse and rat fibroblast cells, as well as several mouse tissue slices for many different stimuli including temperature, PMA, DEX, Forskolin, sodium bicarbonate, LiCl, and CoCl2. This variety of cell types and stimuli provide compelling evidence that the method is broadly useful. This is exciting since it promises to reveal the very insightful PRC with a single experiment from a desynchronized state rather than numerous experiments from a synchronized state for many different phases of the input.
After reading the methods section multiple times and consulting the previous paper on the SR method (Masuda et al. Nat. Commun. 2021), the SR method remains opaque. As I understand it, the method steps are 1) measurements of the individual cell rhythms give estimates of the order parameter R and phase Omega; 2) scale R by beta to account for 'experimental condition' (Omega is left unchanged); 3) assuming the population is desynchronized, invert Eq.(4) to find the PRC g(theta) or more specifically the parameters (a, alpha) in the underlying phase oscillator model that gives rise to the PRC. So that this method can be useful, it is imperative that the exposition of the method be improved greatly. It should be possible for a reasonably competent reader to implement the method for their own data. Here are some specific reasons why I find the method difficult to understand: a) The parameters omega and delta t in Eq.(4) are not defined until much later (after Eq. (12)) in a different context.

Based on your comment, we have added an explanation of the parameters in Eq. (4) in
the section following Eq. (4). (Page 21, b) Possibilities for the PRC g(theta) are given in Eq.(5) and Eq.(6), but those are for fitting to PRC data. The g(theta) resulting from the SR method is implicit in the phase oscillator model Eqs.(9) -(12). Perhaps give the expression for the PRC referenced in (17).
We have added an alternative description of the experimentally obtained PRC g(theta) and g(theta) calculated using the simulation model as exp ( ) and model ( , Δ ) , respectively. c) All of the validation results of the method come before the method is adequately explained in the methods section. Perhaps the introduction could include a better summary of the details of the method or an example using synthetic data (generate 'data' from the phase oscillator model and then use the method to infer the known PRC).

As per your suggestion, we have added a summarized experimental method to the new
Supplementary Figure 1 (Page 3-4, Line 68-96). d) Use the index j instead of i to avoid confusion with the imaginary unit i.
We have corrected the text to avoid this confusion (Page 21, Line 594-600).
e) The SR method is better explained in Masuda et al. Nat. Commun. 2021. For completeness, the details there should be duplicated here.
We have added a description in the revised manuscript (Page 21-23, Line 587-644).
f) The details of the inversion process (step 3 above) are not described. This is a nontrivial root finding problem especially since the integral in Eq.(4) is a smoothing operation.
I am surprised that the authors do not report experiencing any numerical difficulties (and how they were overcome) with this inversion.
We previously showed that when t is small enough compared to the period (at least approximately 1/3 of the period), R monotonically increases with an increase in a, as shown in Supplemental Figure 2a in Masuda et al., Nature Communication, 2021. Thus, we can uniquely determine a using R. We have added a detailed explanation to the supplemental text (Page 8, Line 214-218, Supplementary Fig. 8, Supplementary text).
In short, my primary suggestion is that the clarity of the description of the SR method be improved. The paper would also benefit from more investigation of the mechanism behind and the limits of the SR method. Why does it work? When does it not work?
Thank you for your insightful comments, we have modified the manuscript thoroughly so that our method is easier to understand.
Minor suggestions: 1) Change "mouse" to "mammalian" in the title since the approach was applied to multiple mammalian clocks.
Thank you for your suggestion, we have edited the title.

Other revision:
1) It was difficult to distinguish between a and α in Eq. (10), so we corrected these to a and b.
2) We have corrected the letters of SR amplitude and phase and the SR parameters obtained from PRC to ′ , ′ , ′ and ′ ,respectively.
3) Fig. 1h and Fig. 4  In their reply, the authors have answered all my questions and comments very suitably and convincingly. They have performed complementary experiments which confirmed and also improved the original results. The revised manuscript is improved and more clear. Therefore, I highly recommend this manuscript to be accepted.
Reviewer #2 (Remarks to the Author): The explanation of the SR method is much better in this revision. I especially appreciated the description on lines 587-600 and the new figures in the supplement. The revised manuscript satisfactorily addresses the issues that I previously raised.
minor suggestions: 1) Eq. 4: Should g have a subscript "exp" as appears in the text after Eq. 4? If not, remove "Here, " from the paragraph following Eq. 4.

Responses to reviewer's comments.
We thank all the reviewers for their constructive comments. We have revised the manuscript accordingly.

Reviewer #2 (Remarks to the Author):
The explanation of the SR method is much better in this revision. I especially appreciated the description on lines 587-600 and the new figures in the supplement. The revised manuscript satisfactorily addresses the issues that I previously raised.
minor suggestions: 1) Eq. 4: Should g have a subscript "exp" as appears in the text after Eq. 4? If not, remove "Here, " from the paragraph following Eq. 4.
We have added the subscript "exp" to g in Eq. 4.
2) line 596: "Each PRCs are described as follows," -> "The PRCs have the assumed forms," We have corrected it.
3) Eq. S2: FYI I checked that this formula is correct.
Thank you so much for confirming it. 4) Supp. Fig. 1: "to estimate PRC" -> "to estimate the PRC"; add "the" before PRC and SR in several other places.
We have mentioned this in the figure legend.
We have added a period. Fig. 2de: insert a space between "Concentration" and "(uM)"

7)
We have inserted a space. 8) line 608: "proportional coefficient" -> "The proportionality coefficient beta"; ", and may include" -> " accounts for" We have corrected it. 9) line 217: "the phase of stable point of PRC" -> "the phase of the stable point of the PRC" We have corrected it.