Oscillation detection with period uncertainty and a limited sampling budget

Most rhythm detection experiments involve measurements collected at regular time-intervals. These equispaced designs will fail to detect oscillations at particular acrophases and frequencies even when oscillation amplitude is large. The spurious false negatives pose a challenge for studies aiming to determine the presence or absence of oscillations across a range of parameters. Here, we present a method to construct sampling schemes that are robust to parameter uncertainty, and demonstrate their improvements relative to equispaced designs. We prove that maximizing the worst-case statistical power is equivalent to a mixed-integer quadratically-constrained programming problem with convex relaxation. Using this equivalence, we construct optimal and near-optimal designs for a range of experimental conditions. Our method also allows us to include regularization and sampling constraints, ensuring that improving statistical power via measurement scheduling does not introduce unnecessary complexity in the design.


Introduction
Biological oscillations serve various essential roles in living systems [1].When studying these systems experimentally, measurement timing is crucial for ensuring that oscillations are not overlooked as random noise.In general, protocols which measure at regular time intervals, known as equispaced or uniform designs, are treated as the gold standard for rhythm detection.If the period of the oscillations is known, this preference is justifiable since equispaced measurements have optimal statistical properties for fixed-period harmonic regression [2,3].On the other hand, systems with an unknown periodicity are outside the scope of this justification.In addition, ethical constraints may prohibit the collection of measurements at regular time intervals.For example, a study with human participants may be required to only collect data during waking hours e.g. if the participants are from a vulnerable population.Hence, systems with period uncertainty or sampling constraints pose new challenges in the optimal design of rhythm detection experiments.
When considering period uncertainty, a shortcoming of equispaced measurements can be illustrated in a simulated rhythm detection experiment (Fig. 1).All of the simulated oscillations have the same amplitude, and are randomly assigned an acrophase and a long, medium, or short period.When measurements are collected at regular time intervals, the rhythm detection results are artificially depleted for certain acrophases, even though oscillations are equally present across all acrophases.This depletion occurs because short period oscillators of certain acrophases are always near the same value when measurements are collected (Fig. 1A).On the other hand, if measurements are collected at irregular time-intervals, the short-period oscillations are detected without any evidence of phase depletion, while the medium and long period oscillations continue to be detected reliably (Fig. 1B).
Examining other measurement budgets, phase-dependent power consistently emerges when the frequency of the oscillation is close to half the measurement budget (Fig. 1C).This critical frequency is known as the Nyquist rate of the design.While it is not unusual for spectral features to behave poorly near the Nyquist rate, we investigated whether this loss of performance could be routinely circumvented through irregular sampling.
Using an exact expression for statistical power [4], we investigate how parameter uncertainty and the measurement budget determine the applicability of irregular sampling.We find that the worstcase statistical power, meaning the lowest value of the power across all signals of a given amplitude, is a reasonable figure of merit for this problem.We prove that worst-case power maximization can be reformulated as a mixed integer quadratically constrained (MIQCP) programming problem with convex relaxation.This MIQCP is closely related to Elfving optimality (or E-optimality), a topic in optimal experimental design which has been well-studied for a variety of statistical models including harmonic regression [2,5].
The remainder of this paper is structured as follows.Section 2 provides background necessary for understanding our approach to power maximization.The numerical consequences of our main theoretical result are investigated in Sec.3.1.Using a standard MIQCP-programming method, we perform worst-case power optimization for a range of budgets and period uncertainties.These solutions perform well in worst-case scenarios but introduce unwanted variability in the power.We reduce power variability by adding a regularization term to the objective and enforcing equispaced constraints on a subset of the measurements.In Sec.3.2 we demonstrate that the regularization and equispaced constraints lead to more robust designs by showing that they perform better when periodogram-based hypothesis testing is performed.Sec.3.3 provides a proof of our equivalence theorem.We conclude in Sec.3.4 by investigating convergence properties of the pure-regularization problem.For certain values of the grid coarseness, we find numerical evidence that the irregular candidate solutions are globally optimal.To make this work accessible to all readers, we provide an R package known as Power-CHORD (Power optimization via Convex programming for HOmoscedastic Rhythm Detection) for maximizing power using our equivalence theorem, available from the first author's GitHub [6].Oscillations of all periods can now be detected robustly, even though the measurement budget is the same as in the equispaced case.(C) Even when the measurement budget is large, phase-dependent power emerges near the Nyquist rate when power is plotted as a function of acrophase and frequency.

Power analysis for harmonic regression
Harmonic regression is a popular statistical framework for studying systems with a dominant frequency [7].The single-frequency cosinor model with MESOR (Midline Estimating Statistic of Rhythm) Y 0 , amplitude A, acrophase ϕ, and frequency f is given by in which ε(t) ∼ N (0, σ) is homoscedastic Gaussian white noise.For a known frequency f , Eq. ( 1) can be rewritten as so that all unknown parameters appear linearly.Given data y = {y i } N i=1 measured at distinct times t = {t i } N i=1 ⊂ [0, 1] N , the optimal least squares solution for estimating the coefficients β = (β 0 , β 1 , β 2 ) appearing in Eq. ( 2) The amplitude and acrophase of the signal can be estimated from β, to obtain Â = β2 1 + β2 2 and φ = arctan ( β2 / β1 ).
Using the cosinor model, rhythm detection can be formulated as a hypothesis test in which the null hypothesis is β 1 = 0 = β 2 and the alternative is β 1 ̸ = 0 or β 2 ̸ = 0. Equipped with a formal statistical test, the statistical power can be used to quantify how reliably the experimental design detects oscillations.Since the power of a study is influenced by factors that are unknown prior the experiment, it is useful to construct designs that achieve high power in a range of conditions.
For rhythm detection, assuming a prior distribution on the frequency and a minimum effect size, one would then be able to determine how many samples are necessary to achieve at least a given statistical power.We adopt this framework here, and seek experimental designs that achieve high power for a range of frequencies and acrophases.
Monte Carlo simulation is a simple technique for estimating the power of a proposed study.This type of strategy is too computationally expensive for numerically optimizing the power.Fortunately, the statistical power of the one-frequency cosinor model can be evaluated analytically under the assumption of homoscedastic white noise.
Theorem 2.1 (Power of one-frequency F test, [4]).Consider the one-frequency cosinor model applied to data {y i } N i=1 collected at times (t i ) N i=1 , with amplitude A, noise level σ, frequency f , acrophase ϕ, and type I error rate α ∈ (0, 1).Suppose we test the hypotheses with H 0 as the null and H 1 as the alternative, using the test statistic is the total sum of squares, RSS = ||y − X β|| 2 2 is the residual sum of squares, and r = 3 is the number of degrees of freedom in the model.Then, the power γ(A, σ, f, ϕ) is given by the expression in which is the F-distribution on (n 1 , n 2 ) degrees of freedom, and F λ (x; n 1 , n 2 ) is the noncentral F distribution with (n 1 , n 2 ) degrees of freedom with noncentrality parameter λ ≥ 0. For λ = 0, the noncentral distribution F 0 agrees with the standard F distribution.
Equation ( 3) is significantly cheaper to evaluate than Monte-Carlo-based power estimation.Beyond accelerating our computations, Eq. ( 3) also influences our work due to its mathematical structure.As we will show in Sec.3.3, analysis of Eq. ( 3) reveals a convex optimization problem underlying worst-case power maximization.

Fisher information matrix
The Fisher information matrix (FIM) appears in many areas of statistics, including the theory of optimal experimental design [8].By encoding figures of merit in linear-algebraic properties of the FIM, design problems can be converted to convex programming problems [9].Our work makes use of the FIM in a similar way, and so we recall some basic properties and definitions relevant to our model.

Definition 2.1 (Fisher information matrix). Given a linear regression model
We consider two Fisher information matrices in this paper: the Fisher information matrix corresponding to the one-frequency cosinor model f (t) = [1 cos(2πf t) sin(2πf t)], which we denote by and the Fisher information matrix corresponding to the de-meaned model f (t) = [cos(2πf t) sin(2πf t)], denoted by We refer to M as the reduced Fisher information matrix.When studying the Fisher information matrix in an optimization setting, it is common to restrict the measurement times t to a fixed partition τ = (τ 1 , . . ., τ K ) ⊂ [0, 1] for some large K ≫ N .The definition of the Fisher information matrix can be modified accordingly.
We can use Defn.2.2 to rewrite Eq. ( 5) as in which the weights µ i ∈ {0, 1} obey K i=1 µ i = N .By working with a fixed partition and searching for an optimal choice of weights, the decisions variables are now discrete instead of continuous real numbers.Still, this discrete formulation is preferable because the decision variables µ i appear linearly.

Results
For the remainder of the paper, we scale time so that the duration of the entire study is T = 1 and set the noise strength σ = 1 to reduce the number of parameters in the model.In these units, the amplitude parameter A can be interpreted as a measure of the signal-to-noise ratio.We also follow the convention that the acrophase ϕ ∈ [0, 2π) is expressed in radians so that it is easier to compare phase-dependent features between signals of different frequencies.

Optimized sampling improves worst-case statistical power
We construct two types of irregular measurement schedules: threshold designs and balanced designs.
Threshold designs are constructed by maximizing the worst-case power for a fixed amplitude and a range of frequencies and acrophases.We use the terminology "threshold" because such designs are generated with the aim of scheduling measurements such that the power is above the highest possible threshold for all parameters of interest.In contrast, balanced designs are constructed to achieve a favorable trade-off between worst-case power and variability in power.
Our focus on worst-case power in the construction of threshold designs is motivated by our model of parameter uncertainty.We assume that N measurements are to be collected over the course of an experiment that seeks to investigate all frequencies f in a given window f min ≤ f ≤ f max and no information is available about the acrophase of the signal before the experiment takes place.
Without any knowledge of the acrophase, it is impossible to alter the average power by rescheduling measurements.
Proposition 3.1 (Conservation of average power).For a prior distribution on frequency and acrophase of the form p(f, ϕ) = 1 2π ρ(f ), the average value of the noncentrality parameter is given by in which N is the measurement budget.
Proof.By Theorem 2.1, the timing of measurements only influences the power through the noncentrality parameter given in Eq. ( 4).Integration against a uniform prior shows that the noncentrality parameter is only sensitive to the number of measurements.
Corollary 3.0.1.Provided that the acrophase prior is uniform, the average power with respect to frequency and acrophase depends only on the number of measurements and not when they are scheduled.
For the remainder of this work, we model parameter uncertainty using a uniform prior in phase and frequency Since it is not meaningful to optimize average power in this setting, we focus instead on worst-case power.Our main theoretical result converts worst-case power maximization into an equivalent but better studied optimization problem.We defer the proof of Theorem 3.1 to Sec. 3.3 and focus first on the numerical consequences.
Theorem 3.1.For a rhythm detection experiment with N > 3 distinct measurements and homoscedastic noise, finding a design t * that maximizes power uniformly in acrophase and frequency is equivalent to solving a convex mixed-integer quadratically constrained programming problem in which Q 1 , . . ., Q K are symmetric positive semi-definite matrices.
The optimization problems in Eq. ( 6) and Eq. ( 7) are equivalent in theory, but their structural differences have computational consequences which influence their applicability.Since the cost function in Eq. ( 6) is cheap to evaluate, various stochastic optimization methods [10] would be reasonable choices for constructing near-optimal designs.On the other hand, it may be useful to estimate how close a proposed design is to the theoretical optimal performance.From this perspective, the formulation in Eq. ( 7) is advantageous because it allows the problem to be studied by standard MIQCP programming methods.These methods provide estimates of the difference in performance between a candidate design and the theoretical maximum.Such estimates can be time-consuming to produce, with run-times heavily dependent on the parameters under consideration.We find that after adding additional structure to the optimization problem, we are able to obtain numerical evidence of the optimality of certain candidate designs.Since our work aims to understand the theoretical limitations of power optimization while also constructing high-quality numerical solutions, Eq. ( 7) is best-suited to our purposes.
Threshold designs are constructed by maximizing Eq. ( 7) numerically.When compared in a worst-case scenario, threshold designs can perform significantly better than equispaced designs (Fig. 2).For various measurement budgets, irregular sampling is only worthwhile when the Nyquist rate of the equispaced design is included in the frequency window under investigation.This introduces a threshold in the (f min , f max ) plane, above which one can expect substantial benefits from irregular sampling (Fig. 2A).Examining the distribution of measurement times for a particular experimental setting (N = 48, f min = 1, f max = 24) reveals that measurements are pulled closer together by the optimization method, reducing the temporal coverage of the experiment but improving its ability to detect high frequency oscillations (Fig. 2B).Unlike the equispaced designs which undergo power oscillations at integer multiples of the Nyquist rate, the irregular designs do not have a strong dependence on the acrophase or frequency of the signal (Fig. 2C).At the same time, removal of the strong phase dependence is accompanied by weak power oscillations throughout the acrophase-frequency plane, especially at low measurement budgets.While irregular sampling improves worst-case power, there is a clear need for additional structure in the optimization problem.
We consider two strategies for adding structure to the optimization problem and thereby striking a better balance between worst-case power and power variability.The first strategy involves adding a penalty term to the objective function to introduce a bias away from solutions with a high degree of power variability.This is an example of regularization, a common technique for finding robust solutions to optimization problems.In our problem, we use the average gradient of the noncentrality parameter to construct a regularization term.This term takes the form given in Prop.5.5 and is simplified using Cor.5.0.1.Added to the optimization problem, we are left with an expression of the form in which the weight w reg > 0 determines the strength of the regularization and R is the symmetric matrix given in Eq. ( 26).
The second strategy involves assigning a subset of the measurements to an equispaced measurement grid (Fig. 3).Since equispaced measurements have low worst-case power but appear to exhibit little power variability, we expect that a balance between power in a worst-case scenario and average power variability could be found by searching for designs that contain an equispaced sub-design.
Equispaced constraints also simplify the final design.We use the terminology "balanced designs" to refer to optimal and near optimal solutions to Eq. ( 8) with optional equispaced sampling constraints.
To gauge the influence of regularization and equispaced constraints on balanced designs, we imposed equispaced constraints that require measurements to be collected either every 1 hr, 2 hr, or 3 hr.For a fixed amount of frequency uncertainty (f min = 1,f max = 24), the presence of equispaced constraints, explicit regularization, and the measurement budget together determine the relationship between worst-case power and variability (Fig. 4).As expected, the stricter equispaced constraints force the design to be more similar to an equispaced design and therefore limit performance in a worst-case setting.Explicit regularization slightly diminishes power variability while unconstrained designs generally achieve the best balance between worst case power and variability.Overall, the various balanced designs considered here perform comparably well to threshold designs in terms of worst-case power and power variability.

Periodogram analysis with balanced and threshold designs
Our study of threshold and balanced designs has been limited to the setting of a fixed-period cosinor model.If such designs were to be employed experimentally, the data they produce would likely be analyzed in broader statistical frameworks.For instance, periodogram analysis may be performed as an initial test of periodicity before fitting a fixed-period cosinor model to the data.In this section, we simulate periodogram analysis with threshold, balanced, and equispaced designs to determine if they remain applicable in this more general statistical context.Several approaches to periodogram analysis are commonly employed in the study of periodic data.We make use of the Lomb-Scargle periodogram, a method developed for the study of irregularly sampled data [11,12].Given a list of frequencies, the Lomb-Scargle periodogram returns a power for each frequency which can be interpreted as the intensity of that frequency component in the data.High intensity frequencies, which manifest as peaks in the periodogram, can then be investigated using hypothesis tests.We detect oscillations using a hypothesis test derived from the Lomb-Scargle periodogram [13] which is included in a standard periodogram library [14].
We investigated the reliability of equispaced, balanced, and threshold designs for Lomb-Scargle based rhythm detection study in a simulated dataset.The dataset was made up oscillatory and white noise signals of various amplitudes, acrophases, and frequencies (Fig. 5).Performance differences between designs were most noticeable at the highest and lowest frequencies under consideration.
Importantly, the balanced design performed well at high and low frequencies, whereas the threshold design failed to detect the lowest frequency oscillations and the equispaced design could not resolve the highest frequency oscillations.These differences in performance persisted when the amplitude and abundance of oscillatory signals were increased.A randomly generated design was included for reference and was not reliable at detecting high frequency oscillations, similar to the equispaced design.This supports the argument that the performance improvements in the threshold and balanced designs are stronger than would be expected of a design that deviates randomly from equispaced measurements.

Proof of the equivalence theorem and derivation of the method
We first consider a single frequency and then generalize to multiple frequencies.We aim to choose a set of N measurement times t = (t 1 , . . ., t N ) such that in which γ(t; A, f, ϕ) is the power function given in Eq. ( 3).We first show that the optimization problem in Eq. ( 9) is equivalent to maximizing the noncentrality parameter from Eq. ( 4).Proof.Use Prop.5.1 to express F (x; n 1 , n 2 , λ) as in which z = n1x n1x+n2 .Since n 2 = N − 3 > 0 by assumption, it follows that z ∈ [0, 1] and we may apply Prop.5.2 and Prop.5.3 to justify term-wise differentiation of the series.Differentiate term-wise and apply the same reasoning as in the proof of Prop.5.3 to find which verifies that F (x; n 1 , n 2 , λ) is monotone decreasing for all λ ≥ 0.
It follows from Prop.3.2 that the optimal t * appearing in Eq. ( 9) can be found by simply maximizing the noncentrality parameter.A connection between the non-centrality parameter and the reduced design matrix further simplifies the problem.
Proposition 3.3.Maximizing the noncentrality parameter uniformly over acrophase is equivalent to maximizing the smallest singular value of the reduced design matrix X = [cos(2πf t) sin(2πf t)] .
Proof.The parameter λ can be expressed as We reparameterize Eq. ( 10) to obtain max By homogeneity, we may assume A 2 = 1, reducing our problem to finding The min-max principle for the singular values of a matrix A ∈ C n×2 is Here σ 2 ≤ σ 1 , so i = 2 gives the smallest singular value.When i = 2, the right hand side of Eq. (11) agrees with the square of the right hand side of Eq. ( 12).
To further simplify the optimization problem, notice that M = XT X is symmetric and therefore has real eigenvalues.
Corollary 3.1.1.Since σ 2 2 ( X) = λ min ( M ), maximizing power uniformly in acrophase is equivalent to maximizing the smallest eigenvalue of the reduced Fisher information matrix M .
The remaining steps in reducing our optimization problem are motivated from a computational perspective.Corollary 3.1.1reveals the convex structure of the problem, since maximizing the smallest eigenvalue of a symmetric positive semi-definite matrix is a concave programming problem [15].We simplify further so that the problem can be solved using standard convex programming packages.We also assume for the remainder of the proof that the decision variables have been restricted to a discrete grid, as discussed in Sec.2.2.
Since M is low-dimensional, we can use an exact expression for its smallest eigenvalue and convert eigenvalue maximization to a quadratic programming problem.
if and only if in which M (µ; τ , f ) is the Fisher information matrix associated to the bi-variate linear model and • denotes the Hadamard product.
Proof.The smallest eigenvalue is given by M (µ; τ , f ) can be expressed as and the optimization problem in Eq. ( 13) is reduced to Equation ( 17) is an example of an integer-constrained quadratic programming problem.When considering a range of candidate frequencies (f 1 , f 2 ) simultaneously, Eq. ( 17) generalizes to in which A f is the quadratic form of the previous section corresponding to frequency f .To convert Eq. ( 18) to a form suitable for numerical study, we make use of Prop.5.1 to obtain a linear objective with quadratic constraints.In this final form, the problem can be studied numerically using standard solvers such as [16] or open-source alternatives [17].The same conversion is applicable to the regularized problem.

Numerical evidence for optimality in the regularization-only problem
In this section, we investigate numerical solutions to the regularization-only problem The optimization problem in Eq. ( 19) has fewer quadratic constraints than the problem in Eq. ( 8) because the frequency averaging has been carried out analytically.This difference in structure makes easier to find numerical evidence for the optimality of candidate solutions and affects the degenerate critical points symmetry of the problem (Fig. 6) When solving an MIQCP, it is common to impose a time limit for the computation.In some cases, the solver will terminate early because the value of the objective function at the candidate solution is sufficiently close to the optimal solution to the continuous relaxation.When this occurs, it provides numerical evidence for the optimality of the candidate solution.This type of numerical evidence was occasionally obtained for the balanced designs in the MIQCP studied here, but in most cases, the solver would not terminate until it reached the time limit.
We identified large regions of parameter space where the regularization-only problem halts because the solver has produced evidence of optimality (Fig. 7).The optimal solutions to the problem have a structure that depends on the boundaries of the frequency prior.Dense measurements are collected early in the experiment and depending on the parameters of the prior, sparser measurements may be collected near the end (Fig. 7A).Varying the measurement budget together with the boundaries of the prior reveals regions of parameter space where irregular measurements are optimal (Fig. 7B).For a fixed measurement budget, there is a threshold above which equispaced measurements again appear to be optimal.It should also be noted that the coarseness of the underlying partition is of central importance in producing evidence for optimality.Determining if the evidence of optimality identified here remains valid at finer partitions remains to be understood.Benefits of irregular sampling are dictated by period uncertainty and measurement budget.(A) Threshold designs were generated for uniform frequency priors f ∼ Unif(f min , f max ) and various measurement budgets.Each circle represents the difference in worst-case power between a threshold design constructed by our method and the corresponding equispaced design.Comparing the gain for different values of f min and f max for each budget reveals an uncertainty threshold above which irregular sampling becomes worthwhile.As the budget increases, a greater amount of frequency uncertainty can be tolerated before irregular sampling becomes substantially beneficial.(B) Examples of threshold designs for a high amount of frequency uncertainty (f min = 1,f max = 24).Even at high budgets, the measurements in threshold designs tend to cluster in a sub-portion of the experiment.(C) Influence of frequency and acrophase on power for signals of a fixed amplitude (A = 1).The improvements in power achieved by the threshold design come at the expense of introducing weak power oscillations throughout the acrophase-frequency plane.At a high budget, these effects are negligible suggesting that threshold designs should only be considered when the budget is sufficiently large.13 0 T The remaining measurements are assigned to the unoccupied times in the grid in such a way that the worst-case power is maximized.The result is a design where measurements are to be collected at both the blue and red time points.
Figure 4: Tradeoff between worst-case power and power variability We investigate the tradeoff between worst-case power and power variability assuming a fixed amount of frequency uncertainty (f min = 1, f max = 24) for various measurement budgets and types of optimized design.Constraints mainly affect worst-case power and not the regularization score.The best trade-off is generally achieved by designs that do not adhere to equispaced measurement constraints  1) with acrophases distributed uniformly at random.The parameter p osc dictates the probability that a randomly generated signal in the dataset will be oscillatory.Signals are classified as oscillatory if their periodogram produces a statistically significant p-value after FDR correction.The extent to which oscillations are reliably detected is influenced by their abundance (p osc ), as well as their frequency and amplitude.Detection is summarized visually using an receiver operator characteristic curve for each design.The equispaced design fails to detect high frequency oscillations, whereas the threshold design performs favorably at high frequencies but poorly at low frequencies.The balanced design performs comparably well at all frequencies of interest.For reference, a design points distributed uniformly at random is also included and performs generally worse than the optimized designs.To investigate how discretization affects the optimality, we studied different ratios between measurement budget and discretization (f min = 1).Depending on the budget, discretization, and frequency prior, the MIQCP solver can show that a particular type of sampling is optimal or fail to prove optimality within the provided computation time.Comparing the panels reveals a discretization-dependent region in the frequency-budget plane where irregular sampling is optimal.
Our work has identified conditions where irregular measurement schedules significantly improve the statistical power of rhythm detection in the presence of parameter uncertainty.We provide an efficient numerical method for optimizing the power in such settings by reformulating worstcase power optimization as an MIQCP.We add structure to the optimization problem through regularization and equispaced constraints.These features simplify the designs without significantly worsening their performance, making them more applicable to realistic experiments.Application of our designs to a periodogram-based hypothesis test demonstrates that they remain applicable in statistical frameworks beyond the cosinor model.Focusing on the regularization-only problem, we obtained numerical evidence for optimality of irregular designs at various degrees of frequency uncertainty.
Several immediate extensions of our method could be considered in future work.First, many studies are constrained by windows of time in which measurements cannot be collected.Such constraints could be added to our method without altering its structure.An alternative approach accommodates such constraints using weighted least-squares [18] and could potentially be combined with our framework to allow for simultaneous optimization of measurement times and weights.
Second, statistical power is not the only figure of merit relevant to reliable design of rhythm detection experiments.Depending on the goals of the study, it may be important to understand how the design of the experiment influences the accuracy of period estimation.Various periodogram and spectral techniques are available, many of which have a long history of application to irregularly sampled data [14,19,13,20].Explicitly treating spectral accuracy as an additional figure of merit could further improve usefulness.Third, our method could be applied to adaptive experimental design where the prior distributions are updated using Bayesian inference.Closed form expressions for the regularization kernel would no longer be available, but nothing would be fundamentally changed in making the method adaptive.If the method were extended to other priors, it would be worthwhile to compare the uniform priors considered here to discrete priors made up of Dirac masses centered at integer frequencies.This would still be analytically tractable and could lead to qualitative differences in the optimal solutions of the problem.
The assumption of homoscedastic Gaussian white noise and the cosinor-based hypothesis test are central to our method but also limit its applicability.Extending the method to more general statistical frameworks and understanding its relation to methods in nonlinear optimal experimental design [21] both provide promising directions for future research.Another meaningful extension to the method could come from treating the number of replicates at each time-point as a decision variable in the optimization problem.This is common in optimal experimental design and could provide additional benefits to the statistical power.Through continued investigation in settings including and beyond power analysis, we hope that irregular sampling techniques may join the vast array of tools available to experimental biologists.

Methods
The following propositions are used in Prop.3.2 to verify that the noncentral F-distribution depends monotonically on the noncentrality parameter .
Proposition 5.1 (Series representation of noncentral F-distribution).The cumulative distribution function F (x; n 1 , n 2 , λ) for a noncentral F-distributed random variable with degrees of freedom n 1 , n 2 and noncentrality parameter λ can be expressed as a series Proof.A derivation of this formula is given in [22].Proof.The coefficient B(1, a, b) in Eq. ( 21) can be expressed as dt is the gamma function [23].Use integration by parts and the identity Γ(z + 1) = zΓ(z) to obtain Since z ∈ [0, 1] the second term in Eq. ( 22) is of definite sign, hence To obtain the asymptotic result, rearrange Eq. ( 22) and integrate by parts to obtain Since x z a < 1 for x < z, and for b > 0,z ∈ [0, 1], the following integral is finite we can evaluate the desired limit using the dominated convergence theorem to obtain = 0 for any λ > 0. Uniform convergence of the derivatives can be verified using the Weierstrass M-test.To construct a majorant for T k (z, n 1 , n 2 , λ), notice that ∂ λ T k has critical points at λ = 2k and λ = 0. We know that T k (z, n 1 , n 2 , 0) = 0 and lim λ→∞ T k (z, n 1 , n 2 , λ) = 0 so ).This majorant also leads to a bound on the derivative, since |∂ λ T 0 | ≤ 1 2 T 0 and for k > 0 To prove uniform convergence of k ∂ λ T k , it remains to show that lim k→∞ g k+1 g k < 1.Indeed, we have We know lim k→∞ 1 + 1 k k = e and so lim k→∞ with the last equality justified by the second claim of Prop.5.3.
The following proposition is useful for verifying convexity of the reformulation of power optimization in terms of the de-meaned Fisher information matrix.
Proposition 5.4.For any u, v, w ∈ R n the matrix is symmetric and positive semidefinite.
Proof.From Eq. ( 23) it is clear that C = C T , notice that C can be expressed as a positive combination of rank 1 matrices and thus C is positive semidefinite.
The main goal of our method is to maximize the statistical power across a range of parameters.
This direct approach can produce low-amplitude power variations in the frequency-phase plane.
The next proposition shows how we can use L 2 regularization to remove this unwanted variability without altering the convexity of the problem.
Proof.Since λ(µ) = µ T η(τ ), the gradient can be written as The techniques discussed so far produce an integer-constrained optimization problem with convex relaxation.Solving this problem numerically requires a final simplification, whereby the cost function is simplified at the expense of adding additional constraints to the problem.This approach is an application of "disciplined convex programming" [24], an algorithmic framework for posing convex programming problems in a format suitable for numerical investigation.For our problem, the cost function is simple enough that the reduction can instead be carried out by hand using the techniques below.

Figure 1 :
Figure 1: Power fluctuates near the Nyquist rate of equispaced measurements.(A Left) Examples of signals with long, medium, and short periods relative to the Nyquist rate of an equispaced design with N = 18 measurements.Points represent measurements of the oscillations with noise simulated from a normal distribution.The shortest period oscillations occur at 99% of the Nyquist rate of the design.(A Right) An empirical p-value distribution obtained from a fixedperiod cosinor model (12,600 oscillators are simulated for each period) with the significance threshold α = 0.05 indicated by the dashed line.Oscillations are detected robustly at all phases for the 24 hr and 6 hr period oscillators, but the 2.69 hr oscillations have a p-value distribution heavily dependent on acrophase even though the oscillations are equally prevalent at all acrophases.(B Left) To improve detection of fast oscillations, N = 14 measurements are collected in the first half of the experiment and N = 4 in the second half of the experiment.(B Right) The p-value distributions obtained from the new design no longer suffer from artificial acrophase depletion.Oscillations of all periods can now be detected robustly, even though the measurement budget is the same as in the equispaced case.(C) Even when the measurement budget is large, phase-dependent power emerges near the Nyquist rate when power is plotted as a function of acrophase and frequency.

Figure 3 :
Figure 3: Procedure to construct designs with equispaced constraints (Top) The duration of the experiment is cut into a fine grid indicated by the open circles.The user splits their measurement budget into two groups, indicated here by color.(Middle) Measurements in the red group are assigned to a uniform grid over the entire duration of the experiment.(Bottom)The remaining measurements are assigned to the unoccupied times in the grid in such a way that the worst-case power is maximized.The result is a design where measurements are to be collected at both the blue and red time points.

Figure 5 :
Figure 5: Periodogram analysis with irregular sampling detects simulated oscillations Each panel represents a distinct dataset made up of homoscedastic noise and cosinor signals as in Eq. (1) with acrophases distributed uniformly at random.The parameter p osc dictates the probability that a randomly generated signal in the dataset will be oscillatory.Signals are classified as oscillatory if their periodogram produces a statistically significant p-value after FDR correction.The extent to which oscillations are reliably detected is influenced by their abundance (p osc ), as well as their frequency and amplitude.Detection is summarized visually using an receiver operator characteristic curve for each design.The equispaced design fails to detect high frequency oscillations, whereas the threshold design performs favorably at high frequencies but poorly at low frequencies.The balanced design performs comparably well at all frequencies of interest.For reference, a design points distributed uniformly at random is also included and performs generally worse than the optimized designs.

Figure 6 :
Figure6: Influence of regularization on fixed-phase worst-case power maximization (A) An experiment aims to detect oscillations of known phase ϕ = 0 and unknown frequency using N = 10 equally spaced measurements.Two additional measurements are to be added at times t 1 and t 2 .(Left) If the frequency prior is tight, power is high and the additional two measurements have little influence on the objective.(Middle) As the upper limit of the frequency prior approaches the Nyquist rate, the extra two measurement times have a stronger influence.Importantly, the objective function has several local maxima, illustrating the non-uniqueness.(Right) At a higher level of frequency uncertainty, the objective function oscillates faster and becomes increasingly degenerate as a function of the two additional measurement times.(B) After adding regularization to the objective, there emerges a preference for earlier measurement times, reducing the degeneracy in the problem.

Figure 7 :
Figure 7: Regularization favors irregular measurements.The pure regularization problem corresponds to minimizing the mean squared norm of the frequency-acrophase gradient of the noncentrality parameter.This problem is parameterized by the frequency window under consideration, measurement budget, and spatial resolution of the underlying partition.The number of points in the partition is given by N fine = ⌊cN ⌋ for a constant c > 0. (A) Optimal designs for various choices of f min with all other parameters held fixed (f max = 45, N = 56, N fine = 67).Solutions are similar in structure: measurements are collected as rapidly as possible near the start of the experiment and then, depending on the value of f min , a portion of the budget is allocated to slightly sparser measurements.(B) To investigate how discretization affects the optimality, we studied different ratios between measurement budget and discretization (f min = 1).Depending on the budget, discretization, and frequency prior, the MIQCP solver can show that a particular type of sampling is optimal or fail to prove optimality within the provided computation time.Comparing the panels reveals a discretization-dependent region in the frequency-budget plane where irregular sampling is optimal.