Scalable growth models for time‐series multispectral images

Vegetation indices (VIs) are produced as a combination of different reflectance bands that are captured by multispectral images (MSIs). These indices, such as normalized difference vegetation index (NDVI), are reported to be proxy indicators of photosynthetic activity, plant canopy biomass, and leaf area index. To determine the utility of using VI derived from MSI to model plant growth, random regression (RR) models with linear splines and different orders of Legendre polynomials were applied to data collected (years 2019 and 2020) as part of the Genome‐to‐Fields initiative. Growth curves of maize (Zea mays L.) hybrids were modeled using both NDVI and cumulative NDVI (cNDVI) phenotypes. Due to the difference in MSI recording dates, and sparse overlap in hybrids between years, all the analyses were nested within a year. Results indicate that RR models using Legendre polynomials provide a robust and scalable method for modeling growth curves using phenotypes extracted from MSI; however, RR models using linear splines showed inconsistent convergence. Growth curves estimated using NDVI and cNDVI showed low‐to‐moderate heritability (0.11–0.44) and a range of genetic correlations (−0.15 to 0.97) with grain yield. This study demonstrates the utility of MSI for modeling genetic growth trends, with the best modeling results obtained when using Legendre polynomials and cNDVI.

MSIs repeatedly across the growing season and at different locations will allow breeders to monitor the growth and development of the plants, increasing their understanding of genotype-by-environment interactions on crop performance.
Normalized difference vegetation index (NDVI) is a commonly used dimensionless index that is calculated using infrared and red reflectance bands from an MSI, and it has been used to measure biomass and leaf area index in a number of crops (Hatfield & Prueger, 2010).In addition, recent studies have shown that there exists considerable heritable variation in NDVI when collected at different crop developmental stages (Anche et al., 2020;Krause et al., 2019;Rutkoski, Poland, Mondal, Autrique, et al., 2016).Studies have shown that VIs, such as NDVI, have a high level of association with biomass production during the early stages of heading and grain filling (Hatfield & Prueger, 2010).Furthermore, the seasonality and change of crop vegetation and biomass can be assessed through time-series recording of NDVI values.Several methods have been used to characterize sequential time-series NDVI values such as principal component analysis (Hirosawa et al., 1996).Samson (1993) and Townshend et al. (1985) proposed a skew index and range index to quantitatively analyze discrete time-series NDVI values.The integration of sequential NDVI values through time as a cumulative NDVI (cNDVI) was reported to directly correlate with dry-matter accumulation (Tucker et al., 1981) and has been used as a tool to monitor plant biomass across years.
In recent years, the uses of MSI and image processing platforms have become a common practice in plant breeding (Galli et al., 2020;Herrmann et al., 2020;Lopes et al., 2012).However, robust statistical frameworks for the analysis of repeatedly measured phenotypes derived from MSIs are in the early stages of development.In Campbell et al. (2018), it was shown that random regression (RR) models can be used to improve genomic prediction accuracies for longitudinal complex traits in rice.In animal breeding, RR models are common tools to model repeated measures that are recorded on a continuous scale, such as time or age (Boligon et al., 2012;Brito et al., 2017;Huisman et al., 2002;Kirkpatrick et al., 1990;Lopes et al., 2012;Robbins et al., 2005).As the volume of repeated phenotypes collected through time increases, methods are needed that can effectively model high-dimensional longitudinal data sets.RR models are scalable as the number of time points increases due to the use of covariance functions as opposed to covariance matrices.When Legendre polynomials or linear splines are incorporated, RR can be used to model the variance and covariance of parameters at and among the time points, with the number of parameters estimated scaling with the order of the function rather than the number of time points, resulting in far fewer parameters than models that estimate all pairwise covariances between time points (Kirkpatrick et al., 1990;White et al., 1999;Schaeffer, 2004;Misztal, 2006).

Core Ideas
• Phenotypes from multispectral images can be used to model growth curves.• Normalized difference vegetation index (NDVI) and cumulative NDVI (cNDVI) are heritable and correlated with grain yield.• Random regression (RR) models can be used to model growth curves of crops using phenotypes from multispectral images.• RR models using Legendre polynomials and cNDVI provide robust models of crop growth.
In this study, we aim to investigate approaches for modeling growth curves of maize (Zea mays L.) hybrids using MSI phenotypes and RR models.For that purpose, NDVI and cNDVI were used as MSI phenotypes.RR models using both linear splines and Legendre polynomials were evaluated to identify accurate, robust, and scalable modeling approaches for growth curves based on MSI phenotypes.

Agronomic phenotypic data
The phenotypic data were collected from field trials at Cornell University's Musgrave Research Farm as part of the Genomes to Fields initiative.The trials were planted at research sites in Aurora, NY, in 2019 and 2020.The field trial was planted as a randomized complete block design with two replicates in each year.For 2019 and 2020, plot-level phenotypic data were collected for important end-of-season traits for 612 and 403 hybrids, respectively.Given only 16 hybrids overlapped between the 2 years, separate analyses were run for each year.There were 6 and 12 time points at which MSI data were collected for 2019 and 2020, respectively.Table 1 shows the dates on which the MSI data were collected, the growing degree days (GDDs), and growth stages of the maize hybrid plants for the 2019 and 2020 trials at different developmental stages approximated by the GDDs.A MicaSense RedEdge multispectral camera was mounted on an unoccupied aerial vehicle (UAV) and used to collect images in five reflectance bands: red, blue, green, red edge, and near infrared.The UAV was programed to fly at an altitude of 25-30 m above the ground.From each UAV flight, orthomosaics were constructed using Pix4D (https://www.pix4d.com).Quality control of the orthomosaics was done in three steps: (1) ensuring the reflectance panel was properly registered, (2) checking the overlap statistics of the orthomosaic construction log, and (3) applying nonlocal means denoising.The orthomosaics were then used to calculate summary statistics, such as mean, median, and variance for individual reflectance bands, and VI.NDVI was calculated for each pixel using the following equation: where  NIR is the near-infrared reflectance, and  R is the red reflectance.Mean plot level NDVI was then calculated as the mean of all pixel NDVI values within the plot boundary.Among the different summary statistics that were available for NDVI, mean values for each plot were used in this study.Additionally, the cNDVI at a specific time point, , was calculated for 2019 and 2020 data using the rollmean function of zoo package in R statistical software that takes NDVI values and GDDs at each time point.GDDs were calculated as where  max is the maximum temperature,  min is the minimum temperature, and  base = 10 • C is the base temperature.Any temperature below  base is set to  base before calculating the average.Likewise, the threshold for maximum temperature is usually set at 30˚C, which was never reached in 2019 but reached 33˚C in 2020.The GDDs calculated for each time point were used as time covariates in the RR models.

Genotypic data
Genome-by-sequencing data were available for ∼600 K single-nucleotide polymorphism (SNP) markers for 1557 parental inbred lines.After quality control for missing data and minor allele frequency, ∼122 K SNP markers were used to calculate the additive genomic relationship matrix between parents (VanRaden, 2008), where the additive genomic relationship for the hybrids was calculated by taking half of the parental inbred genomic relationships (Anche et al., 2020).
A single-trait genomic best linear unbiased prediction (ST-GBLUP) model was used to estimate the genetic and residual variances for plot level grain yield, NDVI, and cNDVI at each time point: where  is the vector for the raw phenotype, μ is the overall mean, b is a vector for the fixed effect of replicate,  is the vector of random additive genetic effects of the hybrid to be estimated, X is a design matrix for the fixed effect of replicate, Z is the design matrix for the additive genetic effects, and e is the vector of residuals.It is assumed that,  ∼ (0, σ 2  ), where σ 2  is the additive genetic variance, and G is the genomic relationship matrix between the hybrids (VanRaden, 2008).The residuals were assumed to be distributed as  ∼ (0, σ 2  ), where  is the identity matrix, and σ 2  is the variance of random residual effects.Heritabilities for NDVI and cNDVI values at each growth stage were calculated as a narrow-sense genomic heritability: Genetic and residual correlations of grain yield measures at the end of the season with NDVI and cNDVI values were estimated using multi-trait genomic best linear unbiased prediction (MT-GBLUP) as follows: (5) where y 1 and y 2 are vectors of phenotypic values for traits 1 and 2, μ 1 and μ 2 are overall means,  1 and  2 are vectors of random additive genetic effects,  1 and  2 are the incidence matrices linking 𝑏 1 to 𝑦 1 and 𝑏 2 to 𝑦 2 , 𝐙 1 and 𝐙 2 are the incidence matrices linking 𝑎 1 to 𝑦 1 and 𝑎 2 to 𝑦 2 , and 𝑒 1 and 𝑒 2 are vectors of random residual effects for traits 1 and 2, respectively.It was assumed that where G is the same as in Equation ( 4), and ] ⊗ ).Genetic and residual correlations were calculated using the estimated covariances σ 12 and σ 12 , respectively.

Random regression model
RR models with linear splines (RRLS) and Legendre polynomials (RRLP) were used to fit NDVI and cNDVI values at different time points.Here, the RR models were used to continuously model the (co)variance of NDVI measurements at different time points as a function of time.The random effects at any time point are then calculated as a function of the estimated RR coefficients and standardized measure of GDDs.
The RRLS used in this study is similar to the model used in Misztal (2006).The RRLP model with m order Legendre polynomials can be denoted as φ  (), where  is the standardized time that lies between −1 and 1 (Schaeffer, 2004): where  = 1, … , , where  is the number of time points,  min is the earliest time point, and  max is the last time point.In this study, the GDDs are used as a measure of time.
The RRLP model used to fit NDVI and cNDVI at time point  was where  ∶ is the nth observation on ith hybrid at time t in the jth fixed factor at time t, and   is a fixed effect that accounts for the mean growth curve.For the term ∑   = 0    ∶ , m is the order of Legendre polynomial,   is the mth RR coefficient for the additive genetic effect of hybrid , and  ∶ is the mth order coefficient of Legendre polynomial for ith hybrid at time t in the jth fixed factor.For the term ∑   = 0    ∶ ,   is mth RR coefficient for permanent environmental effect for plot , and  ∶ is a random residual effect.
The RRLS model that was used is described as follows (Anche et al., 2020): where  ∶ is the nth observation on ith hybrid at time t in the jth fixed factor, and   is a fixed effect of replicate by time and accounts for the mean growth curve.For the term ∑   = 1    ∶ ,  is the number of the spline knot points,   are the RR coefficients for the mth knot for the ith hybrid, and  ∶ are the covariables related to time t at knot point m.For the term ∑   = 1    ∶ ,   is the random permanent environmental effect for the mth knot for plot l, and  ∶ is a random residual effect.
In RRLS, linear segments connected at knot points are used to form the regression function (Misztal, 2006;Robbins et al., 2005).This approach is computationally efficient but can be sensitive to the selection of knot points.Among the six time points, the first and the last two time points were used as knot points when RRLS was fit for the 2019 data.For the 2020 data, (the 1st, 3rd, 5th, 8th, 11th, and 12th) were used as knot points among the 12 time points.
In matrix notation, the RR model follows where y is the vector of observations,  is the incidence matrix corresponding to the fixed effect of replicate nested within time,  is the vector of fixed effects,  1 is a matrix of regressors, which are either a function of time using knots for RRLS or Legendre coefficients when RRLP is used,  is a vector of regression coefficients for the random genetic,  2 is the matrix of regressors for each time point, and  is the vector of random permanent environmental effect for each plot.The (co)variance matrix for the random genetic effects () is defined as where σ 2 1 is the genetic variance at the first RR coefficient, σ 2  is the genetic variance at kth RR coefficient, σ 1,2 is the genetic covariance between the first and the second RR coefficients, and G is the same as Equation (4).The (co)variance of the genetic effects at any time points can be calculated by pre-and post-multiplying Σ  by matrices containing the standardized time measurements of interest.The (co)variance matrix for the random permanent environmental effect (pe) is where diagonal elements in Σ  are the variance components for permanent environmental effects for each regression coefficient, the off-diagonals are the covariance of permanent environmental regression coefficients, and I is the identity matrix.As with the genetic effect, the (co)variance of permanent environmental effects at any time points can be calculated by pre-and post-multiplying Σ  by matrices containing the standardized time measurements of interest.A diagonal (co)variance matrix was defined for the random residual effects (e) as follows: where σ 2 1 is the residual variance at time point 1, σ 2 ∶ is the residual variance at time point t, and I is the identity matrix.

Heritability and correlations
Figure 1 shows the heritability estimates for NDVI and cNDVI values at different growth stages in 2019 and 2020.These estimates were obtained using the additive and residual variance estimated from the ST-GBLUP model (Equation 3) and calculated using Equation (4).Heritability estimates for NDVI values for the 2019 data ranged from 0.11 to 0.35, with the lower estimate observed at the last growth stage and the highest at the first growth stage (VT).For cNDVI, heritability ranged from 0.34 to 0.35.For the 2020 data, heritability for NDVI values ranged from 0.12 to 0.44, with the lowest heritability observed at the V9-V11 growth stages and the highest heritability observed at the last growth stage (R6).
For cNDVI, heritability estimates decreased from 0.27 at the earliest observed growth stage to 0.20 at the final time point.Similar results from an ST-GBLUP model run for each year separately were observed in Anche et al. ( 2020), where a considerable amount of heritable variation was observed for NDVI across years.The narrower ranges in heritability for cNDVI are expected given it is a cumulative function of the NDVI measurements at each time point.For the 2020 data, decreasing heritability estimates for cNDVI were observed at the later growth stages, likely due to most hybrids reaching maturity.These results indicate there is heritable variation in NDVI and cNDVI that can be utilized to model genetic effects for growth curves.
In Figure 2, we show the genetic correlation of grain yield with NDVI and cNDVI at different time points in 2019 and 2020.For 2019, the genetic correlations between grain yield and NDVI ranged from moderate (0.39) to high (0.97) at different growth stages, whereas the genetic correlations between grain yield and cNDVI were found to be consistent across growth stages ranging from 0.64 to 0.67.For 2020, the genetic correlation between grain yield and NDVI at different growth stages ranged from −0.15 to 0.81, excluding the time point collected at the V12 growth stage as the MT-GBLUP model failed to converge.Lower and negative genetic correlations were observed between grain yield and NDVI at later growth stages in 2020, which could be because the plants entered senescence at the later growth stages.Genetic correlations between grain yield and cNDVI increased through time ranging from 0.25 to 0.62 in 2020.Several other studies have found moderate-to-high genetic correlations between NDVI and grain yield under different environmental conditions (Hassan et al., 2019;Rutkoski, Poland, Mondal, et al., 2016;Sun et al., 2017b).Our results suggest that NDVI and cNDVI values at different growth stages could be used to improve prediction accuracy for important agronomic traits such as grain yield when used as a secondary trait (Rutkoski, Poland, Mondal, et al., 2016;Sun et al., 2017a;Anche et al., 2020).

Random regression
As mentioned above, RRLS and RRLP functions were fit to model NDVI and cNDVI values at different growth stages using the data from 2019 and 2020.For both 2019 and 2020, the RRLS performed well with no convergence issues when only NDVI was used as a phenotype.However, when cNDVI was used as a phenotype, RRLS did not converge for both 2019 and 2020, likely due to high covariance between knots points.When RRLP was fit, the model converged well for both 2019 and 2020 using either NDVI or cNDVI as a phenotype.Therefore, only results from the RRLP will be presented in this section.
For 2019, RRLP models were fit using both second-and third-order polynomials.Lower order Legendre polynomials were utilized to avoid overfitting as only six time points were collected.The order of the polynomial that best fit the data was selected using the Akaike information criterion (AIC) (Table 2).As shown in Table 2, RRLP fit with a third-order  polynomial provided the best fit (AIC = −4473.86)when NDVI was used as a phenotype, and a second order polynomial provided the best fit (AIC = 6239.15)when cNDVI was used as a phenotype.For 2020, with 12 time points, the RRLP was fit using second-, third-, and fourth-order polynomials.As shown in Table 2, the RRLP with a third-order polynomial had the lowest AIC (= −16793.06)when NDVI was used as a phenotype, and the fourth-order polynomial provided the best fit (AIC = 28630.36)when cNDVI was used as a phenotype.These results suggest that, when the data set contains fewer time points, lower order polyno-mials are robust enough to fit and explain the covariance structure.When the data set contains >12 time points, however, fitting higher order polynomials should be considered.
In addition to the correlations calculated at the recorded time points (Figure 2), NDVI and cNDVI estimates from the RRLP with second-and third-order polynomials were calculated at evenly spaced intervals and correlations with grain yield measured at the end of the season were calculated for each time point using the 2019 and 2020 data sets, respectively.This was done by predicting genetic effects for NDVI and cNDVI at evenly spaced intervals (50 GDD) throughout the growing season using the RRLP.As shown in Figure 3A, the correlations with NDVI were relatively higher at the earlier and later time points and lower at the mid-time points when the 2019 data set was used.On the other hand, higher correlations were observed during earlier relative to later growth stages when cNDVI was used as a phenotype.In 2020, correlations were observed to be lower early and late in the growing season and higher at midseason time points when NDVI was used as a phenotype (Figure 3B).When cNDVI was used as a phenotype, however, the correlation was found to increase through time in the 2020 field season (Figure 3B).This could be because cNDVI accounts for earlier season NDVI data and, therefore, becomes more informative in predicting grain yield toward the end of the growing season.These results are similar to the results presented in Figure 2 for the 2020 data, where correlation between grain yield and NDVI decreased through time, whereas the opposite trend was observed with cNDVI as a phenotype.Although this study only examines two field seasons, our results suggest that using cNDVI as a phenotype to model the covariance structure provides more stable and consistent results compared to using NDVI as a phenotype.
To investigate how well the RRLP models performed, estimates from RRLP models were correlated with the estimates from the ST-GBLUP model where either NDVI or cNDVI values at a given time point were fit in the model.Figure 4 shows the correlation between the genetic estimates from the RRLP model with second (red points) and third (green points) order polynomials and the ST-GBLUP model using cNDVI (A) and NDVI (B) values at different time points from 2019. Figure 4A shows the correlation between the genetic estimates from the RRLP and ST-GBLUP model when cNDVI was used as a phenotype.When cNDVI was used as a phenotype, the correlation between the estimates ranged from 0.98 to 0.99 with a second-order polynomial and from 0.53 to 0.64 when a third-order Legendre polynomial was fit.For NDVI, the correlations ranged from 0.62 to 0.98 (red points) for the second-order polynomial and 0.51 to 0.99 (green points) for a third-order polynomial (Figure 4B).
Figure 5 shows the correlation between the genetic estimates from the RRLP and ST-GBLUP models for the 2020 data set.When cNDVI was used as a phenotype, the correlation between the estimates ranged from 0.96 to 0.99 with a second-order polynomial, from 0.67 to 0.91 with a third-order polynomial, and from 0.40 to 0.96 with a fourth-order polynomial (Figure 5A).When NDVI was used as the phenotype, the correlation between the genetic estimates ranged from −0.05 to 0.97 from the RRLP model with a second-order polynomial (red points), from 0.63 to 0.98 with a third-order polynomial (green points), and from 0.25 to 89 with a fourth-order polynomial (blue points) (Figure 5B).
As observed in Figures 4 and 5, higher correlations between the estimates were observed at all growth stages when a lower order (second-order) polynomial was fit in the RRLP model and when cNDVI was used as a phenotype for the 2019 data set.The same is true for the 2020 data set where higher correlations between the estimates were found when RRLP using a second-order polynomial was fit with cNDVI as a phenotype.These results indicate that, with increasing time points recorded via the MSI, lower order RR models should be robust enough to model the covariance structure between the time points by leveraging all the time points collected.This is especially true for the RRLP model because a convergence issue was encountered with the RRLS model due to high correlations between time points.
Using RR models, growth curves can be calculated for many individuals, partitioning permanent environmental effects from genetic effects.For the 2020 data set, MSIs were collected for more time points compared to the 2019 data set, recording earlier stages of crop growth (V(n)-VT).For that reason, the 2020 data set was used to compare the growth curves estimated using NDVI and cNDVI with grain yield in the hybrids.As shown in Figure 6, hybrids with higher grain yield tend to have higher estimates for NDVI or cNDVI, whereas lower yielding hybrids tend to have lower estimates for NDVI or cNDVI throughout the growth stages.A similar relationship was observed between NDVI and grain yield of maize hybrids in Anche et al. ( 2020).These results indicate that the patterns observed in NDVI and cNDVI throughout the growing season correspond to performance in grain yield.Therefore, growth curves fit using cNDVI or NDVI could be used to understand genetic variation in growth patterns under different environmental conditions, leading to a better understanding of the drivers of genotype-by-environment interactions.
Even with multiple time points recorded throughout the growing period, RRLS had convergence issues with both the 2019 and 2020 data sets.The same was found with MT-GBLUP models where the model often failed to converge when more than two time points were fit in the model.Strong correlations between time points could be one of the reasons for the convergence issues observed with RRLS.The nonlinear increase in parameters makes MT-GBLUP a poor choice for modeling longitudinal MSI data.As the number of time points increase, RRLP proved to be the more robust approach in modeling growth curves and the covariance between correlated time points.

CONCLUSIONS
Scalable modeling of growth curves, both in terms of the ability to collect the phenotypes in a high-throughput manner and analyze data collected at many time points, will enable breeding programs to gain better insights into genetic variation for growth traits and genotype-by-environment interactions.With maize as the test species, we demonstrated that RRLP models are scalable, and growth curves estimated using NDVI and cNDVI are both heritable and correlated with grain yield-an economically important terminal trait.Although these statistical methods provide a robust framework for high-throughput growth modeling, further research is needed to extend these approaches to directly incorporate environmental information and connect growth parameters related to biomass with grain yield.We also acknowledge contributions from Nick Lepak, Josh Budka, Cinta Romay, and other current and past members of the Gore, Buckler, Smith, and Nelson labs.The authors would like to thank the Genomes to Fields consortium for providing phenotypic, genotypic, and image data used in this study.This consortium involves more than 30 researchers representing more than 20 research institutions.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflict of interests.

D A T A AVA I L A B I L I T Y S T A T E M E N T
Details about the initiative and publicly available resources can be found at http://www.Genomes2Fields.org.Agronomic and genotypic data can be downloaded at: https://www.genomes2fields.org/resources/#data-downloads.The summary statistics and plot level vegetative indices can be downloaded at https://imagebreed.org.

F
Heritability estimates for normalized difference vegetation index (NDVI) and cumulative NDVI (cNDVI) at different growth stages in 2019 and 2020.T A B L E 2 Akaike information criterion from the random regression models with different orders of Legendre polynomial for the 2019 and 2020 data sets

F
Genetic correlation of grain yield with normalized difference vegetation index (NDVI) and cumulative NDVI (cNDVI) at different growing stages in 2019 and 2020.

F
Genetic correlation of grain yield with normalized difference vegetation index (NDVI) and cumulative NDVI (cNDVI) in 2019 and 2020.F I G U R E 4 Correlation between genetic estimates from the single-trait genomic best linear unbiased prediction (ST-GBLUP) and random regression (RR) model with a second-or third-order Legendre polynomial when cumulative normalized difference vegetation index (cNDVI) (A) or NDVI (B) was used as a phenotype for the 2019 data set.F I G U R E 5 Correlation between genetic estimates from the single-trait genomic best linear unbiased prediction (ST-GBLUP) and random regression (RR) model with a second-, third-, or fourth-order Legendre polynomial when cumulative normalized difference vegetation index (cNDVI) (A) or NDVI (B) was used as a phenotype for the 2020 data set.F I G U R E 6 Estimated genetic effects of normalized difference vegetation index (NDVI) and cumulative NDVI (cNDVI) for maize hybrids with high, low, or average (mean) grain yield.

T A B L E 1
Planting date, multispectral image (MSI) date, growing degree days (GDDs), and growth stages