Phytoplankton Photophysiology Utilities: A Python Toolbox for the Standardization of Processing Active Chlorophyll-a Fluorescence Data

The uptake and application of single turnover chlorophyll fluorometers to the study of phytoplankton ecosystem status and microbial functions has grown considerably in the last two decades. However, standardization of measurement protocols, processing of fluorescence transients and quality control of derived photosynthetic parameters is still lacking and makes community goals of large global databases of high-quality data unrealistic. We introduce the Python package Phytoplankton Photophysiology Utilities (PPU), an adaptable and open-source interface between Fast Repetition Rate and Fluorescence Induction and Relaxation instruments and python. The PPU package includes a variety of functions for the loading, processing and quality control of single turnover fluorescence transients from many commercially available instruments. PPU provides the user with greater flexibility in the application of the Kolber-Prasil-Falkowski model; tools for plotting, quality control, correcting instrument biases and high-throughput processing with ease; and a greater appreciation for the uncertainties in derived photosynthetic parameters. Using data from three research cruises across different biogeochemical regimes, we provide example applications of PPU to fit raw active chlorophyll-a fluorescence data from three commercial instruments and demonstrate tools which help to reduce uncertainties in the final fitted parameters.

The uptake and application of single turnover chlorophyll fluorometers to the study of phytoplankton ecosystem status and microbial functions has grown considerably in the last two decades. However, standardization of measurement protocols, processing of fluorescence transients and quality control of derived photosynthetic parameters is still lacking and makes community goals of large global databases of high-quality data unrealistic. We introduce the Python package Phytoplankton Photophysiology Utilities (PPU), an adaptable and open-source interface between Fast Repetition Rate and Fluorescence Induction and Relaxation instruments and python. The PPU package includes a variety of functions for the loading, processing and quality control of single turnover fluorescence transients from many commercially available instruments. PPU provides the user with greater flexibility in the application of the Kolber-Prasil-Falkowski model; tools for plotting, quality control, correcting instrument biases and high-throughput processing with ease; and a greater appreciation for the uncertainties in derived photosynthetic parameters. Using data from three research cruises across different biogeochemical regimes, we provide example applications of PPU to fit raw active chlorophyll-a fluorescence data from three commercial instruments and demonstrate tools which help to reduce uncertainties in the final fitted parameters.

INTRODUCTION
The first uses of in-vivo chlorophyll-a fluorescence measurements were as a proxy for chlorophylla concentration of photosynthetic organisms (Lorenzen, 1966). More recent advances in active chlorophyll-a in-vivo fluorescence enable measurements of photosynthetic efficiency (Kolber and Falkowski, 1993). Active chlorophyll-a fluorescence has been widely adopted by the oceanographic community as a rapid, non-destructive technique to assess ecosystem status and microbial function across large temporal (Suggett et al., 2006a) and spatial scales (Behrenfeld et al., 1996;Suzuki et al., 2002;Moore et al., 2005Moore et al., , 2006a approaching those of entire oceanic ecosystems (Behrenfeld and Kolber, 1999;Behrenfeld et al., 2006;Suggett et al., 2006b). A variety of instrumentation types have been utilized by the community, from the early models of pump and probe fluorometers (Falkowski et al., 1986), to the more recent PicoF Lifetime Fluorometry (Lin et al., 2016). One of the most prevalent approaches is fast repetition rate fluorometry (FRRf; Kolber et al., 1998) and its variant fluorescence induction relaxation fluorometry (FIRe; Gorbunov and Falkowski, 2004). This is predominantly due to the versatility of these instruments. Commercial fast repetition rate or fluorescence induction relaxation instruments exhibit large dynamic ranges in detection sensitivity suitable for application to eutrophic and oligotrophic systems (Röttgers, 2007). The ability to measure multiple parameters such as the functional absorption cross section and electron transport kinetics (Kolber et al., 1998;Gorbunov and Falkowski, 2004;Röttgers, 2007); as well as the capacity to collect measurements in the laboratory (Fujiki et al., 2007;Mckew et al., 2013;Schuback et al., 2015), connected to a flow-through aquatic water system (Behrenfeld et al., 2006;Houliez et al., 2017) or in-situ (Moore et al., 2003;Suggett et al., 2006b;Fujiki et al., 2008Fujiki et al., , 2011, including on autonomous platforms (Carvalho et al., 2020), are also advantageous features of this type of instrumentation.
FRRf and FIRe instruments are similar in that they can both initiate Single Turnover Fluorescence (STF) protocols which sequentially close reaction centers in the timescale of one turnover of the first electron acceptor. The resulting measurement from this technique is a characteristic fluorescence yield transient which has a typical rise and decay of fluorescence, termed the saturation and relaxation phases, respectively. It is the modeling of this transient that produces useful parameters that indicate the redox state of electron acceptors and quantifies the transfer of electrons (Kolber et al., 1998). The biophysical model developed by Kolber et al. (1998), known as the Kolber-Prasil-Falkowski (KPF) model, describes the changes in the chlorophylla fluorescence with time recorded in these measurements using the parameters F o , the transient minimum fluorescence, F m , the transient maximum fluorescence, σ PSII , the effective absorption cross-section, and ρ the probability of energy transfer between individual PSIIs, also known as the connectivity coefficient. Furthermore, the model also parameterizes a series of time constants of increasing length (τ 1 , τ 2 , τ 3 , τ PSII ) from the relaxation phase describing electron transport from primary electron acceptor of PSII (Q A ) to photosystem I (PSI) via the secondary electron acceptor (Q B ) and plastoquinone pool (PQ).
The KPF model has been widely adopted for the processing of FRRf measurements by fluorometer users and manufacturers, incorporated into custom-written code (Laney, 2003(Laney, , 2008Barnett, 2007;Ciochetto, 2017;Jesus et al., 2019); and proprietary "out-of-the-box" software (Chelsea Technologies Group, 2012a,b;Photon Systems Instruments, 2019;Kolber, 2021). Whilst these software programs have served the chlorophyll-a fluorescence community well there are significant gaps which need to be addressed in order to achieve a global database of highquality chlorophyll-a fluorescence derived parameters (Lawrenz et al., 2013;Hughes et al., 2018). Noticeable differences exist in the application of FRRf to derive photosynthetic parameters e.g., number of sequences averaged per measurement (e.g., Houliez et al., 2017;Kulk et al., 2018), and application of the KPF model in data processing e.g., single or triple exponential model for relaxation kinetics (e.g., Moore et al., 2006b;Perkins et al., 2018) or choice of fixed vs variable ρ parameter (Suggett et al., 2001;Moore et al., 2003). In addition, different instruments exhibit artifacts in the recorded fluorescence transients that are specific to the instrument make, model and or serial number (Barnett, 2007;Barnett in pers. comm.;Laney and Letelier, 2008). To date, there is no freely available software which can accommodate all of these variations in instrumentation, methodology and data processing.
Over time, the chlorophyll-a fluorescence community has strived to improve best practice to account for known uncertainties in FRRf measurements such as detailed instrument characterization (Barnett, 2007;Barnett in pers. comm.;Laney and Letelier, 2008), blank correction (Cullen and Davis, 2003), spectral corrections for differences in excitation and actinic light wavelengths (Suggett et al., 2001;Silsbe et al., 2015) and optimizing signal-to-noise during the measurement of and fitting of data from low biomass systems (Suggett et al., 2005). However, these practices are still not always adopted in part due to lack of awareness and also the limited tools available for data processing, especially when processing high-throughput datasets. It is common for example for single turnover fluorometer users to report values directly from the software (e.g., Wilson et al., 2015;Houliez et al., 2017). However, very few of the proprietary software programs provide any statistical metrics on the robustness of the fit, leaving users with options to either visually identify robust fits which is unfeasible for highthroughput data, rely on averaging parameters after fitting which can still retain considerable error or ignore the step entirely. Additionally, as new methodology and corrections [e.g., baseline correction of Boatman et al. (2019)] become available and the chlorophyll-a fluorescence community move toward more opensource data, complementary and rapidly adaptable open source software for processing FRRf data will be of significant value to progressing community goals.
In this paper we present an open-source tool, Phytoplankton Photophysiology Utilities (PPU), that can be used to process raw chlorophyll-a fluorescence data from a variety of commercial instruments, where we focused on the Chelsea Technologies Group FAST Tracka I and FAST Ocean and Sea-Bird Scientific (previously Satlantic) FIRe. This paper presents the detailed methods of PPU using example datasets from the three instruments that were deployed in flow-through mode across different biogeochemical regions, demonstrating the versatility and robustness of the optimization routines whilst also providing the necessary statistical metrics to perform quality control. Lastly, we recommend steps toward improving the optimization of chlorophyll-a fluorescence data and highlight future steps to be taken to further enhance PPU for the user community.

Data Collection
Single Turnover Active Fluorescence Data Active chlorophyll-a fluorescence was measured using three instruments across three different cruises (Figure 1) as follows: (1) Chelsea Technologies Group (CTG) FAST Tracka I (SN 05-4845-001): Data were obtained during the Irminger Basin Iron Study (IBIS) on a cruise of the RRS Discovery to the high latitude North Atlantic (D350) from 28 April 2010 to 10 May 2010 ( Figure 1A). The CTG FAST Tracka I was connected to the underway seawater supply to collect measurements in the dark chamber according to the data collection protocol outlined in Table 1. Data was converted from binary to a text file with bin2txt.exe program provided with the Laney V6 code Laney (2003Laney ( , 2008 Figure 1B). The CTG FAST Ocean Instrument was connected to the underway seawater supply to collect measurements according to the data collection protocol in Table 1 using a combination of LEDs; blue (λ450 nm), blue+green (λ450+530 nm) and blue+orange (λ450+624 nm). Data was exported to a CSV file using the FastPro8 copy all function and pasted into an excel file. For the purposes of this analysis only data from the blue LED (λ450 nm) is reported.  Figure 1C). The FIRe instrument was connected to the underway seawater supply to collect measurements according to the data collection protocol outlined in Table 1, the protocol consisted of only using the blue LED (λ455 nm) for the excitation flashlets. Original raw data files were exported with no conversions applied.

Particulate Absorption Data
Samples for particulate absorption were also collected and used to demonstrate spectral correction features of the Phytoplankton Photophysiology Utilities toolbox (see section "LED Spectral Correction"). Following the quantitative filter pad technique (Mitchell et al., 2002), 0.5-2 L of seawater were filtered onto 25 mm GF/F, snap frozen in liquid nitrogen and stored at −80 • C for later analysis in a Shimadzu UV-2501 spectrophotometer. Optical density (OD) was measured using the internally mounted integrating sphere approach using an air baseline (Stramski et al., 2015) from 350 to 750 nm (1 nm resolution). Blank filter measurements were also collected. Phytoplankton particulate absorption coefficients were calculated with the Phytoplankton Photophysiology Utilities toolbox as explained in section "LED Spectral Correction."

PPU Analytic Software and Optimization Algorithms
A custom analysis package, Phytoplankton Photophysiology Utilities (PPU) was developed to robustly analyse data from CTG FAST Tracka I, CTG FAST Ocean and Sea-Bird Scientific (previously Satlantic) FIRe instruments. More instruments are supported but will not be discussed further here. The package is available open source on the GitLab repository and free for use following the MIT license 1 ; the package can also be downloaded from the python package index (pip) with the command: pip install phyto_photo_utils. PPU was written in Python (V3.7; 2 ) using robust optimization methods from the SciPy package (available at 3 ). Additional dependencies include tqdm (available at 4 ), numpy (Harris et al., 2020), pandas (available at 5 ), matplotlib (available at 6 ) and sklearn (available at 7 ). The default method in the PPU toolbox for least squares optimization uses a similar approach to Laney (2003), applying a trust region reflective algorithm Li, 1994, 1996) which reduces the potential of the estimated 1 https://gitlab.com/tjryankeogh/phytophotoutils 2 www.python.org 3 www.scipy.org 4 www.tqdm.github.io/ 5 www.pandas.pydata.org/ 6 www.matplotlib.org/ 7 www.scikit-learn.org/stable/index.html parameters becoming "trapped" in local minima or at the user defined boundaries. Further defined parameters for the loss function, tolerance limits and number of function evaluations before the termination allows the user to balance computational efficiency against computational time. The optimization routines fit the Kolber-Prasil-Falkowski biophysical model (KPF;Kolber et al., 1998) with a series of user defined functions and options for the saturation phase and the relaxation phase of the fluorescence induction curve (details below).

PPU Statistical Parameters for Data Quality Control
The PPU fitting routines return a suite of statistical parameters which assess the accuracy, bias and precision of the optimization routines and derived parameters, some of which, but not all are included in the output of manufacturer and custom written software for processing FRRf and FIRe data. Other software programs report the coefficient of determination (R 2 ), χ 2 and reduced χ 2 as goodness of fit tests for the robustness of the non-linear model (Supplementary Table 1). However, we recognize that these metrics are not suitable tests for non-linear models (Andrae et al., 2010;Spiess and Neumeyer, 2010) or one-hit Poisson functions similar to the KPF model (Laney, 2003) and have chosen to exclude these metrics. The following statistical metrics were calculated for the model performance: root mean squared error (RMSE, Equation 1, Chai and Draxler, 2014) assessing accuracy between the fitted model and observed fluorescence yield; normalized RMSE where the mean squared error is divided by the mean fluorescence values (nRMSE, Equation 2); and normalized bias (modified from equation 3 of Seegers et al., 2018), calculated as Equation 3, indicating the systematic direction of the average difference or error between the model and observations. The precision of the fitted parameters is also assessed through fit errors, one standard deviation error for each derived parameter computed with an inverse jacobian covariance matrix normalized to the mean squared residuals. In cases of F o and F m the fit error values are returned as percentages so as to avoid any specific instrument bias on the parameter error interpretation. In each equation O i and M i are the nth flashlet of the observed and modeled fluorescence yield, respectively, and Ô is the mean observed fluorescence.

Functions of the Package
There are 17 functions available within the PPU toolbox (Figure 2), which are described in Supplementary Table S2.
Detailed explanations of each function are given in the documentation online 8 with examples in the GitLab repository demo Jupyter notebook 9 . Here we will demonstrate the application of the primary functions and assess their performance against select proprietary manufacturers software distributed with the instruments and custom-written open-source software from the chlorophyll-a fluorescence community.

Saturation Phase Fitting
The PPU toolbox offers three separate models to fit the saturation curve within the saturation.fit_saturation function, varying in the application of the connectivity coefficient ρ. Model 1, "no ρ" model (Equation 4; Equation 1 in Laney, 2010), assumes no connectivity between photosystems. Model 2, "ρ" model (Equation 5 and 6; Equation 1 from Kolber et al., 1998) assumes connectivity between photosystems and iteratively obtains a value for ρ. Model 3, "fixed ρ" model (Equations 5 and 6) also assumes connectivity between photosystems, but ρ is fixed with a user defined value ("fixed ρ" model). The "fixed ρ" model is commonly applied to datasets collected in low biomass ecosystems, where the chances of accurately estimating ρ are decreased (Suggett et al., 2001), see below for more details. The PPU toolkit includes fixed, and user set bounds to constrain the saturation parameter values. Fixed constraints are set for F o as the intercept of a linear fit of the first eight fluorescence values ± 10%, for F m as the intercept of the of the last 24 fluorescence values ± 10%, and ρ (where applicable) as 0 ≤ ρ ≤ 1. The functional absorption cross section σ PSII ( • A 2 RCII −1 ) bounds are user set (see Supplementary Table S3 for default bounds).

Relaxation Phase Fitting
To parameterize the time constants of PSII reoxidation kinetics from the "relaxation" phase of a fluorescence induction curve, the PPU toolbox relaxation.fit_relaxation function can apply either a single exponential equation (Equation 7; Laney, 2003Laney, , 2010 or a triple exponential (Equation 8; Kolber et al., 1998) to parameterize the decay in fluorescence. The relaxation parameters are also constrained by bounds, with F oRelax confined to the mean of the last three fluorescence values ± 10%, and F mRelax between the mean of the first three fluorescence values ± 10%. In instances where the fluorescence decay of the first few flashlets of the relaxation phase is very rapid, a user defined number of flashlets from the end of the saturation phase of the fluorescence induction curve can be included to the mean calculations determining F mRelax bounds to improve the estimate of F mRelax . The relaxation time constants τ 1 , τ 2 , τ 3 and τ PSII (thought to be equivalent to electron transport from , and corresponding amplitudes (α 1 , α 2 , α 3 ) are constrained by user defined values in µs and default values are presented in Supplementary Table S3.

Blank Correction
Fluorescent materials within the water column, such as colored dissolved organic material and other pigmented detrital material, although not capable of variable fluorescence, can still artificially inflate the background fluorescence yield resulting in an increase in F o and F m , and hence F v /F m (Cullen and Davis, 2003). The PPU toolbox can be used to calculate the blank fluorescence value from raw data files, under the assumption that background fluorescence has no variable yield. A blank correction can then be applied either during the fitting of the fluorescence induction curve or post-fitting by directly correcting F o and F m and F v /F m using the PPU calculated value or user set values. As the requirement for blank corrections has long been considered a best practice requirement (Cullen and Davis, 2003), the effects of blank correction will not be discussed further.

LED Spectral Correction
The spectra of the light emitting diodes (LED) supplying the excitation and actinic irradiances within the FRRf and FIRe instruments do not directly represent the natural available light field for phytoplankton photosynthesis in the water column or the action spectrum of photosynthesis for phytoplankton (Emerson and Lewis, 1942;Neori et al., 1986;Szabó et al., 2014). Failure to correct for these spectral differences can lead to over/underestimations of σ PSII (Suggett et al., 2001;Szabó et al., 2014) and hence spectral correction is now considered best practice, where possible (Hughes et al., 2018). The PPU toolbox employs spectral correction using the most common approximation of the spectral dependency of light absorption by photosystem II (Suggett et al., 2004;Silsbe et al., 2015;Schuback et al., 2017), the phytoplankton specific particulate absorption spectrum [a phy (λ)]. Other approximations include the use of a phy (λ) corrected for photosynthetic vs non-photosynthetic pigments and the fluorescence excitation spectrum derived from passive or active spectral fluorescence measurements (Silsbe et al., 2015). PPU provides two functions to achieve this goal; calculating the phytoplankton specific particulate absorption coefficients from measurements of optical density and calculating spectral correction factor (SCF) for correction of σ PSII and the actinic light (where used) such as in the CTG FAST Act system used for fluorescence light curves.
Following the IOCCG best practice protocols for spectrophotometric measurements of particulate absorption using filter pads (IOCCG Protocol Series, 2018), the PPU function calculate_chl_specific absorption provides options to subtract blank filter measurements and calculate total particulate absorption from the optical density as per Equation 9 (Equation 5.3 and 5.4 in IOCCG Protocol Series, 2018) where OD f is the blank corrected optical density of the sample filter (absorbance; with option to normalize to infrared wavelengths), V is the volume of sample filtered (m 3 ), A is the effective (sample) area of the filter (m 2 ) and β 1 and β 2 are the coefficients for PPU: Processing Chlorophyll-a Fluorescence Data from the LaneyV6 code (y-axis) versus PPU (x-axis) from the D350 FAST Tracka I dataset from the North Atlantic. The "ρ" model was used in both fitting routines. Frontiers in Marine Science | www.frontiersin.org the pathlength amplification accounting for scattering through the filter (see Stramski et al., 2015 for the method appropriate coefficients to apply).
If users have depigmented absorbance measurements (i.e., optical density of non-algal particles after removal of pigments with an extraction solvent) then these measurements can be subtracted to yield phytoplankton specific particulate absorption a phy (λ). Alternatively, the non-algal contribution to the total particulate absorption measurements can be algebraically calculated using an iterative best-fit approach of Bricaud and Stramski (1990) to determine the absorption by nonalgal particles and hence a phy (λ). If phycobiliprotein pigment concentrations or other accessory pigment concentrations such as chlorophyll-b or chlorophyll-c concentrations are high (or suspected) in the measured sample, the selection of the green wavelength absorption (at 580 nm) may not satisfy the criteria that absorption by pigments must be minimal at 580 nm so that the ratio with 692 nm is close to 1. The PPU offered solution is to find the green wavelength between 580 and 600 nm which best satisfies the requirements, although it is acknowledged that this is not a perfect solution. Finally, the calculated a phy (λ) can be normalized to a user set chlorophyll-a concentration to yield the chlorophyll-specific phytoplankton absorption coefficients a * phy (λ)(m 2 mg chla −1 ).
The function calculate_instrument_led_correction calculates the spectral correction factor for correcting either σ PSII or the actinic light as needed by the user. The function utilizes the phytoplankton particulate absorption coefficient spectra [a phy (λ), can also be chlorophyll specific a * phy (λ)] calculated as above or input by the user. For correcting σ PSII as Equation 10, the instrument LED spectra [E LED (λ) 400-700 nm] and background light E background (λ, 400-700 nm) as either the instrument actinic light E actinic (λ, 400-700 nm) or in-situ light field E in−situ (λ, 400-700 nm) are required. For correcting the actinic light both E actinic (λ,400-700 nm) and E in−situ (λ,400-700 nm) are required (Equation 10). The E LED (λ, 400 − 700 nm) for the FIRe, FAST Trackka and FAST Ocean instruments and E actinic (λ, 400-700 nm) for the FAST Act system used in this study are provided in example files with the PPU package. As the measurement of in-situ light fields (400-700 nm) are not always common practice in field studies the PPU toolbox can calculate a representative in-situ light field (spectral distribution of irradiance 400-700 nm) using a modified approach of Stomp et al. (2007b,a) and Schuback et al. (2015). The in-situ light spectrum 400-700 nm at depth z is estimated as Equation 12 using a typical incident solar spectrum (400-700 nm from Stomp et al., 2007b,a) and the spectrally dependent light attenuation coefficient [k d (λ)]. Following Morel et al. (2007), k d (λ)(Equation 13a and 13b) is equal to attenuation by pure seawater (a w (λ) from Pope and Fry (1997) and 1 2 b w (λ) from Morel (1974)) and attenuation by phytoplankton, non-algal particles and cultured dissolved organic material as k bio calculated as a function of chlorophyll-a using the coefficient χ(λ) and exponent c(λ). The spectral correction factors for σ PSII or the actinic background light are then calculated as Equations 10 and 11, respectively, where E LED (λ), E in−situ (λ), E background (λ), E actinic (λ) and a phy (λ) have been normalized to the maximum value of their respective spectra:

Optimization in Low Biomass Systems
In low biomass waters the signal-to-noise ratio can significantly impact the accuracy of the optimization routine fitting fluorescence induction curves and the robustness of the fitted results (Suggett et al., 2005). This is of particular importance when operating the instruments in an underway mode and sampling across vast biogeochemical regimes with large dynamic ranges in biomass concentrations. To improve the signal-tonoise, convergence of the model fitting and accuracy of results, fluorescence transients can be merged or averaged before fitting using the tools.remove_outlier_from_time_average function.
As there are large discrepancies in sampling protocols across the active chlorophyll-a fluorescence community, this function employs merging of fluorescence induction transients within a user set time window rather than achieving a set number of acquisitions per measurement average. An additional option is to exclude transients within the time window that exceed a user defined number of standard deviations of the mean.

FIRe Measurement Bias Correction
Benchtop FIRe instruments use a reference excitation profile measured with fluorescent dye in place of direct measures of the incident LED excitation intensity to normalize the fluorescence emission and generate the fluorescence yield (Satlantic LP, 2010). Manufacturers of the FIRe instrument, Satlantic (now Sea-Bird Scientific; note that the instrument is now discontinued) generally supplied two instrument-specific factory default reference excitation profiles for use with different gain settings. However, noticeable artifacts appear in FIRe fluorescence yield where there is a mismatch between the reference excitation profile and actual incident excitation irradiance. Artifacts in the saturation phase may include a large difference in between the fluorescence yield of the first and second flashlet (Supplementary Figure 1A) as the LEDs are warming up after being turned on, which in some extreme cases can be more than 100% greater than the average flashlet difference. As the excitation protocol transitions from the saturation to relaxation phase, a second significant mismatch occurs between the reference profile and actual incident LED irradiance. This results in a large difference between the last flashlet of the saturation phase and the first flashlet of the relaxation phase (Supplementary Figure 1A) possibly due to a time lag during the protocol switch, as well as a significant drop in fluorescence yield of the relaxation phase (far below F o of the saturation phase) roughly equal to the difference between the last flashlet of the saturation phase and the first flashlet of the relaxation phase (Z. Kolber in pers. comm.). These biases have significant impact upon the fitted results which will be discussed later (see Section "Measurement Bias Correction"). In the release of custom software code FIReworx (Barnett, 2007) for processing FIRe data, Barnett (now Ciochetto) advocated strongly for individual users to collect their own repeated reference excitation profile measurements with fluorescent dye Rhodamine Blue to account for the mismatch during the transition from saturation to relaxation phases of the measurements and to drop the first flashlet (or first few flashlets) from the saturation fitting process. In the case that custom reference excitation profiles are not available or where custom profiles still do not correct the bias the PPU toolbox provides options to correct the bias. Within the saturation.fit_saturation function users can utilize options to ignore flashlets from the start of the saturation phase during fitting. In addition to improving estimates of F mrelax via the relaxation.fit_relaxation function as detailed in section "Relaxation Phase Fitting, " users can employ tools.correct_fire_instrument_bias to correct the fluorescence yield of the relaxation phase by adding the difference between the first flashlet of the relaxation phase and a user set position in the relaxation fluorescence yield (e.g., second flashlet) to all flashlets in the relaxation phase (except the first flashlet; Supplementary Figure 1B).

Demonstration of PPU Functions
To demonstrate the sensitivity, reliability and advantages of the PPU toolbox for fitting single turnover fluorescence transients, data from the FIRe, FAST Tracka I and FAST Ocean were fit using the PPU toolbox and available manufacturers software or custom open-source code.
To demonstrate saturation phase fitting by the PPU toolbox, the "ρ" saturation model was selected within the saturation.fit_saturation function to process data from the FIRe, FAST Tracka I and Fast Ocean . FIRe and Fast Ocean data processing from PPU were compared with manufacturer software, FIRePro (V1.3.2) and FastPro8 (V1.0.55, supplied with the benchtop FIRe and Fast Ocean instruments, respectively). As the analytic software originally supplied with the FAST Tracka I (FRS.EXE) is no longer supported by the manufacturer, data from this instrument was processed with the publicly available custom fluorescence yield analysis software, V6, developed for MATLAB by Laney (2003) and an iteration of the software provided (in pers. comm.) by C.M. Moore (presented in Moore et al. (2007)) which applies the same fitting routines but also reports fit errors. To demonstrate relaxation phase fitting by the PPU toolbox, FAST Ocean data was processed with the relaxation.fit_relaxation function using both the triple exponential decay model to derive τ 1 , τ 2 and τ 3 , and single decay model to derive τ PSII . Both models were applied in order to compare to FastPro8 (v1.0.55) as the software only returns a single decay value τ which is akin to τ 1 . All data from the FIRe and FAST Ocean were blank corrected prior to fitting. Blank measurements were not available for the FAST Tracka I data. The upper and lower bounds for the various derived parameters were not changed from the preset in the Laney V6 code and modified Laney code (see Supplementary Table S4 for bounds and initial estimates) and could not be changed within proprietary software FastPro8 (v1.0.55) and FIRePro (V1.3.2). After fitting by PPU or other software/code, modest filtering was performed to remove data points where F v /F m was > 0.65 or F v /F m < 0, ρ < 0.001 or ρ > 1 and σ PSII < 0 or σ PSII > 2,000 and any of the returned τ parameters were equal to the set upper or lower limits.

Statistical Comparisons Between PPU and Other Software
Parameters derived from the PPU toolbox and various manufacturers and custom software were contrasted and assessed using linear least-squares regression, the mean absolute relative difference [Software MARD, Equation 14 (Gerbi et al., 2016)] and the Software Model Bias (Equation 15, Seegers et al., 2018) where M1 i and M2 i are the nth derived parameters from two different software/code fitting routines. Two variants of the least-squares optimization algorithm were also compared for performance, the trust region reflective and Levenberg-Marquardt algorithms using the Software MARD and Software Model Bias statistics. The Software MARD provides an assessment of the overall error (or difference) between the outputs from two different softwares, whereas the Software Model Bias provides a simple description of the direction of the error between the two software outputs.

(14)
Software Model Bias = 10 Additional functions within the PPU toolbox are presented here to highlight best efforts to reduce uncertainties in processing single turnover fluorescence data. The derived value of the connectivity coefficient ρ can have a significant impact upon other retrieved parameters. To investigate the effect large uncertainties and hence the value of the connectivity coefficient ρ on F v /F m and σ PSII , a series of sensitivity tests were performed where the connectivity coefficient was fixed at set values (0.1-0.9) using the saturation.fit_saturation function with "fixed ρ" model and the results were compared to the recommended value of ρ = 0.3 (Suggett et al., 2001). The correction of biases in the FIRe dataset and attempts to improve the signal to noise is demonstrated through application of the saturation.fit_saturation option to ignore the first flashlet and the tools.remove_outlier_from_time_average function, modifying the set time window from 3-30 min and excluding values ± 3 S.D of the mean within the time intervals, and fitting the fluorescence transients with the saturation.fit_saturation "ρ, " and the resulting change in RMSE and fit errors are reported. Finally spectral correction of σ PSII is demonstrated on the same FIRe dataset using the _spectral_correction.calculate_chl_specific absorption and _spectral_correction.calculate_instrument_led_correction functions.

Laney Code Versus Phytoplankton
Photophysiology Utilities (FAST Tracka I, North Atlantic, D350) Retrieved parameters, F v /F m and σ PSII from the PPU toolbox and Laney (2003) V6 code (referred to as LC from here, Figure 3) for the FAST Tracka I dataset from the North Atlantic (D350) were highly linearly correlated (R = 0.99 and 0.96, respectively), with slopes close to 1 for F v /F m (0.94; Figure 3A) and σ PSII (1.18; Figure 3B). The PPU code provides more statistical assessments than the V6 code, including the root mean squared error (RMSE), normalized RMSE, the bias of the fit as well as the fit errors for the retrieved parameters. The retrieved results had mean and median RMSE values of 0.006 and 0.004 for D350, with a mean bias of 0.05. Fit errors averaged 0.005 ± 0.004 (2.6% mean error) for F o , 0.0007 ± 0.0004 (0.2% mean error) for F m , 21.27 ± 13.00 (5.1% mean error) for σ PSII and 0.13 ± 0.10 (33% mean error) for ρ. Average fit errors from the modified LC (MLC) were similar 0.006 ± 0.004 (3% error) for F o , 0.001 ± 0.001 (0.4% error) for F m , 19.9 ± 12 (4.4% error) for σ PSII and 0.11 ± 0.06 (38% error) for ρ and the modified LC results were also highly linearly correlated with PPU results (Supplementary Figure 2). The F v /F m and σ PSII Model Bias analyses demonstrate no distinctive biases (9% underestimation to 12% overestimations: Figures 4D,E) between PPU and LC, MLC processing codes. Compared to PPU, the LC processing code slightly overestimated F v /F m by 5% and σ PSII by 12%. The Software Model MARD of 3-13 % for F v /F m and σ PSII across all algorithms and codes (Figures 4A,B) reflects the scatter around the 1:1 lines observed in Figure 3. However, for ρ the Software Model Bias ranged from a 55% underestimation to a 40% overestimation in PPU vs other processing codes ( Figure 4F) and the relative error (Software MARD) far exceeded 100 % (Figure 4C).
Overall PPU produced similar estimates of F v /F m and σ PSII from the FASTTracka I (North Atlantic, D350) dataset as compared to the Laney Code. There are substantial differences in the estimates of ρ, which is explored further in the discussion. The mean and median RMSE of FAST Ocean saturation fluorescence transients from the South Pacific (IN2016_T01) fit with PPU were 0.0064 and 0.0058, respectively. Mean fit errors were reasonable for F o (0.004 ± 0.002, 2.3 % mean error) and F m (0.001 ± 0.003, 0.75 % mean error), σ PSII (48.67 ± 52.46, 13% mean error). However, the fit errors for ρ were considerably higher (0.25 ± 0.15, 269 % mean error).
Generally, there is good agreement between derived parameters from FastPro8 and PPU (n = 4,244 σ PSII , n = 4,250 for F v /F m ) with correlation coefficients for least-squares linear regression of 0.98 and 0.90 for F v /F m ( Figure 5A) and σ PSII (Figure 5B) respectively, with slopes ranging from 0.81 to 0.96. Derived estimates of ρ were poorly correlated between FastPro8 and PPU (Figure 5C; n = 3473). The lower sample size for ρ is due to the FastPro8 not always returning a value. The Software MARD between F v /F m from PPU and FastPro8 was 9% with FastPro8 typically 3% higher (Figure 5D). For σ PSII the Software MARD was 4% and FastPro8 produced values 11% lower ( Figure 5D). There was a strong positive model bias toward higher ρ values reported by PPU (7.01; Figure 5E), i.e., FastPro8 values were 50% lower. Moreover, the values of ρ ranged from 0.002 to 0.371 within the FastPro8 results, whereas they ranged from 0 to 0.99 within PPU. As a result, the Software MARD between PPU and FastPro8 was > 100%.
The fitting of the relaxation phase to derive the decay time constant τ (µs) was performed using both the single and triple decay method as the FastPro8 software only returns a single decay τ value for the initial part of the relaxation phase (akin to τ 1 ). Both the single decay τ PSII ( Figure 6A) and triple decay τ 1 (Figure 6B) values from PPU were compared to the FastPro8 τ 1 . The decay parameters derived from PPU were poorly correlated (Figures 6C,D), although FastPro8 consistently produced larger values of τ 1 compared to the triple decay τ 1 values from PPU. The triple decay τ 1 values are also significantly different from the FastPro8 τ 1 values (Wilcoxon Signed-Rank test, p > 0.05). The single τ PSII values from PPU spanned a wider range of 110-49,046 µs than the FastPro8 τ 1 values of 157-1,840 µs ( Figure 6A). It should be noted that the sample size of the FastPro8 results is much lower (n = 133) than the PPU (n = 2,171) output, suggestive of an issue with the manufacturer software to derive τ 1 .
The PPU processing of FastOcean data from the South Pacific (IN2016_T01) also produced similar values of Fv/Fm and σPSII as compared to the FastPro8 software, but estimates of ρ were again very different, likely due to ρ apparently being constrained to an upper limit of 0.4 in the FastPro8 software. Estimates of the decay time constant τ also varied substantially between PPU and FastPro8, due to the differences in the fitting equation, and is discussed further in the discussion.

Uncertainties in the Estimation of the Connectivity Coefficient (Fast Ocean , South Pacific, IN2016_T01)
The value of the connectivity coefficient, ρ, in FAST Ocean data from the South Pacific (IN2016_T01), demonstrated a significant impact on the value of other retrieved parameters from the saturation phase of fluorescence transients, specifically F v /F m (Figures 7A,C,E) and σ PSII (Figures 7B,D,F). As compared to results when ρ = 0.3, percentage differences reached up to ∼55% for F v /F m and up to ∼35% for σ PSII when the value of ρ increased. This resulted in a Software MARD that ranged from 4 to 56% for F v /F m ( Figure 7C) and 1 to 6 % for σ PSII (Figure 7D). Software Model Bias for σ PSII was ∼1 (Figure 7E) but for F v /F m it decreased with increasing ρ to a minimum value , and ρ (C) from the FastPro8 software (y-axis) versus PPU (x-axis) processing of FAST Ocean data from the South Pacific (IN2016_T01). Bar charts of (D) software MARD (%) and (E) software model bias calculated between FastPro8 and PPU. Note that the Software MARD for ρ in panel (D) is out of the axis range. of 0.72 (Figure 7F), representing a 28% underestimation relative to ρ = 0.3. Note that here, the "software" compared is PPU where ρ is fixed at values 0.1-0.9 as compared to PPU where ρ = 0.3. Uncertainties related to the estimation of ρ are explored further in the discussion.

FIRePro Versus Phytoplankton Photophysiology Utilities
Derived parameters F v /F m and σ PSII from the FIRe dataset from the Southern Ocean (ACE) processed by PPU do not agree well with the results from the FIRePro software (Figure 8). Linear correlations between the PPU and FIRePro derived parameters were poor at 0.12 and 0.28 for F v /F m ( Figure 8A) and σ PSII (Figure 8B) respectively. Note that no further optimization for signal-to-noise or instrument bias have been applied via either PPU or FIRePro. The results show that both retrieved parameters are much higher from the FIRePro software versus PPU, with many of the F v /F m results from FIRePro being above the theoretical maximum of 0.65 for single turnover protocols (Kolber and Falkowski, 1993). Fit errors reported in the PPU processed data were extraordinarily high at 37,805 ± 366,998 for F o , 259 ± 9,970 for F m , 17,241 ± 96,285 for σ PSII and 137 ± 695 for ρ. Overall estimates of F v /F m and σ PSII did not agree well between PPU and FIRePro.

Measurement Bias Correction
The PPU option to correct a bias in the fluorescence transients of the FIRe dataset from the Southern Ocean (ACE) was applied, removing the first flashlet from the saturation fitting (see Supplementary Figure 1). The measurement bias correction decreased the median RMSE from 24 to 19 for the saturation phase fits (Figure 9A). This correction resulted in a shift in F v /F m values (Figure 9B), decreasing the mean values from 0.25 ± 0.08 to 0.21 ± 0.07 (∼16%), whilst also decreasing the mean σ PSII values ( Figure 9C) from 272.9 ± 235.8 to 219.8 ± 180.9 (∼19%). Moreover, there was a decrease in the fit errors, with a 72% decrease in σ PSII error, an 84% decrease in F o error and a 28% decrease in F m error.

Pre-optimization Transient Averaging
There was a change in the density distribution of the RMSE of fitted transients' values as a result of increasing the time window over which the raw transient data is averaged or merged (Figure 10A), which resulted in a decrease in both the mean and median RMSE values ( Figure 10B) and a decrease in the interquartile range (IQR; Figure 10B). For example, an increase in the time average to 10 min decreased the IQR by 133% decrease, which also led to a 53% reduction in the mean and median RMSE values. This shift in the relative robustness of the fit estimated from the RMSE values comes at a cost of the sampling resolution though, an 81% decrease in the number of measurements. The result of this is that a post-processing FIGURE 6 | Density histograms of τ (µs) from the fitted results of FastPro8 τ 1 and PPU using the single decay τ PSII (A,C) and triple decay (τ 1 ) (B,D) fitting routines using FAST Ocean data from the South Pacific (IN2016_T01). Correlation plots of FastPro8 τ 1 (Y-axis) and PPU (X-axes) using the single decay (C) and triple decay (D) fitting routine.
average will result in F v /F m values that are 6% lower than preprocessed averages (Supplementary Figure 3A) and σ PSII values that are 16% lower (Supplementary Figure 3B). A threshold was identified however, at 21 min, after which wider time windows result in increasing RMSE and IQR values.

Spectral Correction
The spectral correction factor (SCF) from the ACE voyage (Southern Ocean; FIRe) was calculated for 197 stations (Figure 11A), averaging 0.60 ± 0.04 (range = 0.46-0.77). To determine its effect upon a σ PSII , an average of a ± 20-min time window was created for each station time point, this average σ PSII,455 was then multiplied by the SCF to create σ PSII, in−situ (Figure 11B). Across the time series there was an average difference between σ PSII,455 and σ PSII, in−situ of 90.83 ± 65.55 ( • A 2 RCII −1 ) ( range = 32.54-464.92; Figure 11C).

DISCUSSION
The aquatic single-turnover chlorophyll-a fluorescence user community has evolved considerably over the last two decades thanks to the growing commercial availability of user-friendly Fast Repetition Rate (FRRf) and Fluorescence Induction and Relaxation (FIRe) fluorometers. Whilst an array of proprietary software and open-source processing codes have served the FRRf/FIRe user community well for deriving photosynthetic parameters from fluorescence transients, the goal of global databases of high-quality single turnover fluorescence data (Lawrenz et al., 2013;Hughes et al., 2018) is not quite achievable with existing data processing tools. In particular no one software program in the current suite available provides adequate statistical metrics for quality control of data, allows for flexible application of the Kolber-Prasil-Falkowski (KPF) model, provides capacity for high-throughput processing (i.e., from continuously recorded datasets) and/or processes data from multiple models/brands of single turnover fluorometers. We have created PhotoPhytoUtils (PPU) to fill these critical gaps in single turnover fluorescence data processing and quality control in a bid to assist the chlorophyll-a fluorescence community to achieve its goals and progress our understanding of the natural variability of aquatic photosynthesis. PhytoPhotoUtils (PPU) produced similar estimates of saturation phase photosynthetic parameters as compared to the FastPro8 software and open-source code V6 from Laney (2003) and modified version from Moore (presented in Moore et al. (2007)) used to process the FAST Ocean and FAST Tracka I data, respectively. The greatest difference was in the estimation of the connectivity parameter ρ which unlike PPU, was not constrained in the V6 and modified V6 code (See Supplementary Table 4) FIGURE 7 | Connectivity coefficient (ρ) sensitivity tests performed on the Fast Ocean data from the South Pacific (IN2016_T01). The percentage difference in F v /F m (A) an σ PSII (B), Software MARD of F v /F m (C), Software MARD of σ PSII (D), Software Model Bias in F v /F m (E) and Software Model Bias in σ PSII (F) from ρ fixed at values 0.1-0.9 as compared to ρ = 0.3. Note that here, the "software" compared is PPU where ρ is fixed at values 0.10.9 as compared to PPU where ρ = 0.3.   The interquartile range (IQR), median and mean RMSE from the different time windows. All data presented here has been corrected for the FIRe instrument saturation bias, i.e., skipping the first flashlet. Data is from the Southern Ocean (FIRe, ACE). and the upper and lower bounds set for the algorithm solution are very likely to impact the final value of ρ. The difficulties in deriving accurate estimates of ρ have been well documented elsewhere (Suggett et al., 2001(Suggett et al., , 2004Bruyant et al., 2005;Babin, 2008), especially in waters with low chlorophyll biomass. The ranges of chlorophyll concentrations in the surface during each campaign presented in this study were 0.5-5 µg/L in the Irminger Basin (Ryan-Keogh et al., 2013), <0.1-0.3 µg/L in the South-West pacific and East Australian Current (data not published) and 0.01-3.7 µg/L in the Southern Ocean (Antoine et al., 2019). Little is known by the authors about the fitting processes and bounds set for parameters within the FIRePro software which produced very different final parameters from FIRe data as compared to PPU output. The software and the desktop FIRe are no longer supported by the current owner company (Sea-Bird Scientific). The fit errors reported by PPU clearly demonstrate that without quality control, or correction beyond removing the blank influence, that the errors for all derived parameters are concerningly high and so it is not surprising that derived parameters from PPU and FIRePRO are poorly correlated. The decay time constants τ 1 , and τ PSII reported by PPU and FastPro8 were also poorly correlated, however, it is acknowledged that the fitting bounds set for these parameters, which are unknown for FastPro8, may strongly influence the final reported values complicating any comparison.
Allowing for the flexible application of the KPF model and modeling of the decay phase of single turnover fluorescence transients is important so that derived parameters can more accurately be compared to other processed datasets, and also to understand the natural variability of and uncertainties in the connectivity parameter ρ and the decay kinetics (τ). Open-source processing codes V6 (Laney, 2003), modified V6 (Moore et al., 2007) and fireworx (Barnett, 2007) provide these options for fitting FAST Tracka (V6 and modified V6) and FIRe (fireworx) PPU: Processing Chlorophyll-a Fluorescence Data data. Our analysis of the uncertainties in the derived value of ρ and subsequent impact of the value of F v /F m and σ PSII supports results from previous studies such as Babin (2008) that showed varying ρ between 0 and 0.6 can lead to variations in F v /F m of ∼20%. As compared to parameters derived when ρ was fixed at 0.3, the value of F v /F m varied up to 55 % and σ PSII up to 35% when ρ varied between 0.1 and 0.9. As suggested by Suggett et al. (2001Suggett et al. ( , 2004) applying the KPF model to saturation phase transient fits with ρ at a fixed value of 0.3 or constrained between 0.1 and 0.4 reduces the impact on F v /F m and σ PSII . The value of τ, the time constant of PSII reoxidation, in nature is poorly understood. A quick Web of Science and Google Scholar literature search found 171 papers using single-turnover fluorometry to study photosynthetic parameters of phytoplankton between 1998 and 2019 and only 40 papers reported values of τ, with variations in the use of single exponential model and triple exponential models to parameterize the decay in fluorescence and the actual time constant reported (i.e., τ 1 , τ 2 , τ 3 , and τ PSII ). The upper and lower limits set for τ 1 , τ 2 , τ 3 , τ PSII in PPU (See Supplementary  Table S3) were constrained using these few examples. Based on the dataset collected from the South-west Pacific Ocean and East Australian current (FAST Ocean dataset), the value of τ PSII fell comfortably within the prescribed bounds of 100 µs to 50 ms (<5% of returned values were equal to the upper or lower bounds). It would appear though that τ 2 is much more rapid (64% of returned values were the lower limit of 800 µs) and τ 3 much longer (70% of returned values were the upper limit of 50 ms) than previously reported. The PPU package will provide the opportunity for the community to explore this variability in the time constants for the reoxidation PSII across a wide range of instruments and hence datasets.
Phytoplankton Photophysiology Utilities has incorporated a number of functions to improve the accuracy of singleturnover fluorescence transient fitting and to reduce biases and uncertainties resulting from background non-variable fluorescence ("blanks"), instrument specific characteristics and high noise in data from low chlorophyll biomass systems. Although not always employed, the correction for nonvariable fluorescence from background dissolved fluorescent material ("blanks") and spectral correction of σ PSII are generally encouraged as standard practice (Cullen and Davis, 2003;Suggett et al., 2004), and so not discussed further here. Other open-source processing codes have drawn attention to the impact of various instrument specific artifacts in driving uncertainties in photosynthetic parameters and have provided relevant corrections (Laney, 2003(Laney, , 2008Barnett, 2007). As we have demonstrated here, failure to correct for the mismatch in the reference excitation profile within the desktop FIRe instrument will result in inflated F v /F m and σ PSII values and high uncertainty in the estimates. When collecting observations in low-chlorophyll biomass ecosystems, protocols to maximize signal-to-noise such as repeated measurements or the averaging/merging of fluorescence transients before fitting must be employed. It has been suggested that uncertainties in derived parameters when signal-to-noise is low (or the absorption cross section is small) could be reduced through averaging of derived parameters after fitting (Suggett et al., 2005;Oxborough, 2013). However, using the desktop FIRe dataset as an example (See Supplementary Figure 3) there is a significant reduction in both the estimated values of σ PSII and uncertainty of the estimates when transients within a 10-minute window are merged before fitting versus the averaging of the derived parameters after fitting. It must be recognized though that even after applying the various corrections to the fitting of single-turnover fluorescence transients that are available in PPU, considerable uncertainty still remains in the estimate of photosynthetic parameters, which is largely ignored by the single-turnover fluorescence community often because such information is not readily available. The mean fit errors of saturation phase parameters F o , F m , and σ PSII derived from PPU ranged from reasonable values of 1-33% in the FAST Tracka I and FAST Ocean datasets to 31-1,050% percent in the FIRe dataset after bias correction and 9-217% percent after bias correction and transient merging within a 10-min interval. With the additional statistical metrics such as fit errors of the derived parameters and root mean square error of the fitting process, users of PPU can rapidly identify and remove data points with high uncertainty, significantly improving the quality of the processed dataset. When selecting a time window with which to average the user must remain aware of the different light histories that phytoplankton will be experiencing, i.e., day and night, and the length of time within a sampling system. Phytoplankton Photophysiology Utilities is an open-source interface between python and active single-turnover chlorophylla fluorescence data that provides automated and standardized data processing for the chlorophyll-a fluorescence community. The PPU package has been developed with the intention that it evolves with community needs and knowledge that will improve the quality and usefulness of chlorophyll-a fluorescence data. As an open-source package, it is easily adaptable to incorporate recently published and future solutions and corrections, for example the "baseline" correction (Boatman et al., 2019). Planned future updates to the package will include support for the loading and processing of data from other fluorometers e.g., Satlantic in-situ FIRe, Satlantic miniFIRe, PSI Fluorometer and CTG LabSTAF, and the export of processed parameters to climate forecasting compliant netCDF format. Some additional functions currently exist in the package which were not featured here, including the fitting of photosynthesis vs irradiance models to fluorescence light curve data. As addressed above, the natural variability in reoxidation kinetics of the PSII-PSI electron chain is not well explored and consultation with the chlorophyll-a fluorescence community will be needed to improve the PPU relaxation kinetics functions and constrain the uncertainties in τ parameters. The PPU package provides increased capacity for users to take greater responsibility for the quality of their single-turnover fluorescence data and achieve community goals of building global databases of high-quality single turnover fluorescence data.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
TR-K and CR conceptualized the idea, wrote the python code, tested the code, collected the data for ACE, and wrote and edited the manuscript. TR-K collected the data for D350. CR collected the data for INV2016_T01 transit voyage. Both authors contributed to the article and approved the submitted version.

FUNDING
This work was undertaken and supported through CSIR's Southern Ocean Carbon and Climate Observatory (SOCCO) programme (http://socco.org.za). This work was supported by CSIR's Parliamentary Grant (SNA2011112600001) and supported by the Australian Government through the Australian Research Council's Discovery Projects funding scheme (project DP DP160103387).