HUB: a method to model and extract the distribution of ice nucleation temperatures from drop-freezing experiments

The heterogeneous nucleation of ice is an important atmospheric process facilitated by a wide range of aerosols. Drop-freezing experiments are key for the determination of the ice nucleation activity of biotic and abiotic ice nucleators (INs). The results of these experiments are reported as the fraction of frozen droplets f_ice(T) as a function of decreasing temperature and the corresponding cumulative freezing spectra N_m(T) computed using Gabor Vali's methodology. The differential freezing spectrum n_m(T) is an approximant to the underlying distribution of heterogeneous ice nucleation temperatures P_u(T) that represents the characteristic freezing temperatures of all INs in the sample. However, N_m(T) can be noisy, resulting in a differential form n_m(T) that is challenging to interpret. Furthermore, there is no rigorous statistical analysis of how many droplets and dilutions are needed to obtain a well-converged n_m(T) that represents the underlying distribution P_u(T). Here, we present the HUB (heterogeneous underlying-based) method and associated Python codes that model (HUB-forward code) and interpret (HUB-backward code) the results of drop-freezing experiments. HUB-forward predicts f_ice(T) and N_m(T) from a proposed distribution P_u(T) of IN temperatures, allowing its users to test hypotheses regarding the role of subpopulations of nuclei in freezing spectra and providing a guide for a more efficient collection of freezing data. HUB-backward uses a stochastic optimization method to compute n_m(T) from either N_m(T) or f_ice(T). The differential spectrum computed with HUB-backward is an analytical function that can be used to reveal and characterize the underlying number of IN subpopulations of complex biological samples (e.g., ice-nucleating bacteria, fungi, pollen) and to quantify the dependence of these subpopulations on environmental variables. By delivering a way to compute the differential spectrum from drop-freezing data, and vice versa, the HUB-forward and HUB-backward codes provide a hub to connect experiments and interpretative physical quantities that can be analyzed with kinetic models and nucleation theory.

1 Introduction

Ice nucleators (INs) of biological and abiotic origins present in aerosols are responsible for facilitating the heterogeneous freezing of atmospheric water droplets above the homogeneous nucleation temperature (Murray et al., 2012; DeMott et al., 2016, 2003). The potential of these aerosols as ice nuclei has significant implications for cloud properties and precipitation patterns (Gettelman et al., 2012; Mülmenstädt et al., 2015; Froyd et al., 2022). Freezing experiments are key sources of information to determine the range of temperatures over which INs promote ice nucleation. The most common method to characterize INs is through immersion freezing experiments, for which a wide range of assays and instruments have been developed. A comprehensive report of various drop-freezing techniques can be found in Miller et al. (2021). The assays are typically performed by placing uniformly sized water droplets with a known IN concentration or area on a substrate or in a multiwall plate that is gradually cooled from a temperature above 0 ^∘C until all droplets are frozen (Kunert et al., 2018; Budke and Koop, 2015). Droplet freezing is detected visually or through the measurement of the latent heat release (Stratmann et al., 2004; Budke and Koop, 2015; Kunert et al., 2018; Reicher et al., 2018), allowing the assignment of a heterogeneous nucleation temperature to each droplet. Drop-freezing experiments record the fraction of frozen droplets, f_ice(T), as a function of decreasing temperature; for soluble or dispersible INs f_ice(T) curves are typically collected at various 10-fold dilutions of the IN sample.

Historically, there have been two interpretations of the dispersion of nucleation temperatures in heterogeneous freezing experiments. The first approach suggests that the stochastic nature of the nucleation process dominates the variability in freezing temperatures (Bigg, 1953; Carte, 1956), while the second approach assumes that the dispersion in temperatures mostly arises from a distribution of nucleation sites (Fletcher, 1969), each with a deterministic, singular nucleation temperature (Levine, 1950; Vali and Stansbury, 1966). Variability in the temperature, volume, and amount of ice-nucleating particles per droplet can also contribute to the dispersion of freezing temperatures (Vali, 2019; Knopf et al., 2020). There is consensus now that both stochastic effects and sample heterogeneities contribute to the distribution of freezing temperatures, and both approaches are used for the modeling of drop-freezing experiments (Vali, 1971; Marcolli et al., 2007; Niedermeier et al., 2011; Murray et al., 2011; Broadley et al., 2012; Wright and Petters, 2013; Herbert et al., 2014; Harrison et al., 2016; Alpert and Knopf, 2016; Vali, 2019; Fahy et al., 2022b). Stochastic modeling of the freezing curves is based on predicting the survival probability of liquid water containing INs as a function of supercooling, and it requires a model for the temperature dependence of the nucleation rate of the IN components. These models have been solved numerically or evolved with Monte Carlo simulations to interpret or resolve the distribution of ice nucleation properties of minerals (Marcolli et al., 2007; Murray et al., 2011; Broadley et al., 2012; Wright and Petters, 2013; Herbert et al., 2014; Harrison et al., 2016) and organics (Zobrist et al., 2007; Alpert and Knopf, 2016) and to perform parametric bootstrapping of experimental data (Wright and Petters, 2013; Harrison et al., 2016). The advantage of the stochastic modeling approach is that it enables a direct link to microscopic properties of the nuclei and can account for the cooling rate dependence of the f_ice(T) data. These approaches require the use of analytical models for the freezing rates and their distribution in the sample.

The modeling of freezing experiments based on the singular approach is based on the framework proposed by Vali (1971). He assumed that each particular IN has a characteristic ice nucleation temperature that is independent of the cooling history. This implies that the IN with the highest characteristic nucleation temperature in a droplet is responsible for its freezing. Given a total number of droplets N₀, the number of frozen droplets N_F(T) at a temperature T gives the range of characteristic freezing temperatures that determines the ice nucleation activity and is used to produce the cumulative freezing spectrum (Vali, 1971, 2014, 2019):

$$\begin{array}{}\text{(1a)}& N_{\mathrm{m}}\left(T\right)={\displaystyle \frac{\mathrm{1}}{X}}[\ln N_{\mathrm{0}}-\ln N_L(T\left)\right]=-{\displaystyle \frac{\mathrm{1}}{X}}\ln [\mathrm{1}-f_{\mathrm{ice}}\left(T\right)],\end{array}$$

where $N_L\left(T\right)=N_{\mathrm{0}}-N_{\mathrm{F}}\left(T\right)$ is the number of unfrozen droplets; $f_{\mathrm{ice}}\left(T\right)=N_{\mathrm{F}}\left(T\right)/N_{\mathrm{0}}$ is the fraction of frozen droplets at temperature T; and X is a normalization factor per unit volume of water, unit mass, or surface of the INs (Vali, 2019). For soluble INs, the normalization factor is commonly defined by the mass of the ice-nucleating material $X=\mathit{\rho }\left(V_{\mathrm{drop}}/\phantom{V_{\mathrm{drop}}d}d\right)$, where ρ is the density of the initial solution, V_drop is the droplet volume, and d is the dilution factor (Kunert et al., 2018). The IN surface area per drop, X=A_drop, is sometimes used as the normalization factor for insoluble INs (e.g., dust, crystals), resulting in a cumulative spectrum per area denoted as N_s(T). However, it is challenging to measure the total IN surface area accurately (Knopf et al., 2020). We note that Eq. (1a) can be used even when the absolute concentrations or areas of the INs are unknown, provided that the user knows the relative concentration of the dilution series derived from a parent sample. The differential freezing spectrum n_m(T) is obtained by the differentiation of the cumulative spectrum (Vali, 1971):

$$\begin{array}{}\text{(1b)}& n_{\mathrm{m}}\left(T\right)={\displaystyle \frac{\mathrm{d}{N}_{\mathrm{m}}\left(T\right)}{\mathrm{d}T}}=-{\displaystyle \frac{\mathrm{1}}{X{N}_L\left(T\right)}}{\displaystyle \frac{\mathrm{d}{N}_L\left(T\right)}{\mathrm{d}T}}.\end{array}$$

The differential spectrum identifies the density of IN active at each temperature and was identified by Vali as the central quantity that can be derived and interpreted from drop-freezing experiments (Vali, 1971, 2019).

The determination of the differential spectrum from the cumulative one by finite differentiation is subject to significant noise, requiring a careful selection of the temperature intervals and extensive sampling (Vali, 2019). As stochastic effects are not considered in the singular temperature formalism, the cumulative and differential spectra should – in principle – depend on the cooling rate (Vali, 1994). The stochastic nature of ice nucleation, combined with the uncertainties associated with the experimental measurements (e.g., different droplet volumes, inhomogeneous samples, different detection efficiencies), can produce significant variations in the cumulative freezing spectra that result in large uncertainties in n_m(T) and we provide its associated Python code and user manual (https://github.com/Molinero-Group/underlying-distribution, last access: 5 May 2023), published in de Almeida Ribeiro et al., 2023). Parametric and nonparametric bootstrapping based on the singular approximation and Monte Carlo simulations have been used to estimate confidence intervals in freezing spectrum measurements (Vali, 2019; Fahy et al., 2022a, b).

A central assumption of the singular freezing approximation is that the freezing of a droplet containing multiple INs is promoted by the IN with the highest nucleation temperature (Levine, 1950). The extreme value sampling is apparent in the concentration dependence of f_ice(T) in experiments (Marcolli et al., 2007; Budke and Koop, 2015; Kunert et al., 2018; Lukas et al., 2022). Using probability theory, Joseph Levine demonstrated that if the distribution of ice nucleation temperatures of the IN population follows an exponential distribution, then the sampling of droplet freezing temperatures corresponds to a Gumbel distribution, and the median freezing temperature T_MED of the droplets scales with the logarithm of the number (or total nucleating area) of IN per droplet (Levine, 1950). Richard Sear more recently demonstrated that Levine's approach is a particular solution for a generalized extreme value problem and used modern extreme value statistics to derive the scaling of T_MED with the number of IN sites per droplet for the three generalized extreme value (GEV) distributions: Gumbel that would arise from an underlying IN distribution with exponential tails, Frechet from those with power law tails, and Weibull from those with an upper cutoff in the freezing temperature of the INs (Sear, 2013). However, there are limitations for the use of the analytical approaches of Sear and Levine for the interpretation of actual drop-freezing data. First, the extreme value sampling results in one of the three GEV distributions only in the limit of an extremely large number of INs per droplet, while in experiments the sampling is typically performed over dilutions down to a few INs per droplet. There is no analytical formulation for the dependence of the extreme value distribution in the low to intermediate concentration regime. Second, the analytical theory assumes that the sampling is complete (i.e., the number of droplets is extremely large), while experiments are typically performed with tens to hundreds of droplets. Third, Sear notes that there is no general analytical theory to predict the GEV distributions from a mixture of populations of nuclei with different temperature dependences (Sear, 2013). In this study we overcome these three limitations through a numerical implementation of extreme value statistics for the modeling of drop-freezing experiments.

A consequence of extreme value sampling is that the differential spectrum n_m(T) represents the underlying distribution of ice nucleation temperatures of all INs in the sample, which we denote as P_u(T), only when the sampling of INs in the drop-freezing experiments is complete. The underlying distribution P_u(T) is akin to a hub that connects the experimental freezing temperatures to physical analysis based on nucleation theory or kinetic and equilibrium models that can elucidate the mechanisms and origins of the distributions of INs (Fig. 1). We here call the cumulative spectrum N_m(T) obtained through Eq. (1a) in this complete sampling limit the intrinsic cumulative spectrum of the system, I_u(T) (Fig. 1). While there is consensus that the quality of the freezing spectrum increases with the number of droplets, a rigorous analysis of how many droplets and IN dilutions should be measured to provide accurate freezing spectra is still lacking. The first goal of the present study is to provide a strategy to optimize the sampling of drop-freezing experiments to derive interpretable differential spectra that are a good approximant of the underlying distribution of heterogeneous ice nucleation temperatures of the sample.

Figure 1 Diagram illustrating the usage of the HUB code: n_m(T) is obtained from the sparsely sampled f_ice(T) or N_m(T) through HUB-backward, and the effect on f_ice(T) or N_m(T) is obtained from the complete sampling of the underlying distribution P_u(T) through HUB-forward. The intrinsic cumulative spectrum I_u(T) is proportional to ${\int }_{T_{\mathrm{m}}}^T{P}_{\mathrm{u}}\left(T^{\prime }\right)\mathrm{d}{T}^{\prime }$ (Sect. 2.2).

The existence of subpopulations or classes in the population of INs (e.g., different classes of bacterial INs, different ice-nucleating sites on complex materials like dust) (Turner et al., 1990) is common in atmospheric aerosols. While several studies have broadly defined populations from the cumulative spectra by the range of nucleation temperatures they encompass (Turner et al., 1990; Creamean et al., 2019) or the origin of the sample (Steinke et al., 2020), there is currently no simple procedure to identify and quantify subpopulations or classes from cumulative freezing spectra N_m(T). The second aim of our study is to map the cumulative freezing spectrum N_m(T) into the differential spectrum n_m(T) in terms of subpopulations that may correspond to different physical nucleation sites in the sample.

To reach the aims above, we develop a method we name HUB (for heterogeneous underlying-based) to model and interpret the results of drop-freezing experiments and provide its associated Python code and user manual (https://github.com/Molinero-Group/underlying-distribution, last access: 5 May 2023). Our method relies on the singular interpretation of freezing experiments: we assume that each individual IN has a characteristic nucleation temperature independent of its cooling history and that the freezing of a droplet containing multiple INs is promoted by the IN with the highest nucleation temperature. This second assumption allows the use of extreme value statistics (Castillo et al., 2005; David and Nagaraja, 2004; Gumbel, 2012; De Haan and Ferreira, 2006) to model and interpret the data.

We present two implementations of the HUB analysis code. The HUB-forward code allows the user to postulate an underlying distribution of heterogeneous nucleation temperatures P_u(T) in the system of interest. The HUB-forward code uses the singular approximation and extreme value statistics to generate an artificial IN dilution series similar to those obtained in experiments, from which it computes the fraction of frozen droplets f_ice(T) and from these derive N_m(T) using Vali's equation (Fig. 1). The HUB-backward code works in reverse, extracting the differential spectrum n_m(T) from a given cumulative N_m(T) using a stochastic optimization procedure (Fig. 1). HUB-backward allows the decomposition of the total population from n_m(T) into subpopulations. The combination of HUB-forward and HUB-backward allows for an analysis of the sensitivity of N_m(T) to the number of droplets and dilutions, as well as the impact of the sampling on the closeness of the differential spectrum n_m(T) to the underlying distribution P_u(T). The determination of distributions obtained from the HUB-backward code could further enable the interpretation of the experimental ice nucleation spectra with the size and structure of INs using nucleation theory, kinetic models, and molecular simulations. For example, Schwidetzky et al. (2023) illustrate the use of the distribution of freezing temperatures obtained with HUB-backward, together with classical nucleation theory for finite surfaces, to interpret the size of the INs of Fusarium acuminatum.

This paper is organized as follows: Sect. 2 presents the methodology, and Sect. 2.1 discusses the details on the implementation of HUB-forward, while Sect. 2.2 describes the HUB-backward procedure to find the differential spectrum n_m(T) and discusses how to determine whether or not n_m(T) has converged to the underlying distribution P_u(T). Section 3 presents examples of applications of both HUB-forward and HUB-backward codes and their capabilities. Section 3.1 analyses the effect of the number of droplets sampled on the cumulative freezing spectrum N_m(T). Section 3.2 uses HUB-backward to compute the differential spectra n_m(T) of various biological INs with increasing grades of complexity in their cumulative freezing spectra. Section 3.3 demonstrates how to extract n_m(T) from the experimental fraction of ice f_ice(T) and the impact of the cooling rate on n_m(T). We end in Sect. 4 with a discussion of the main conclusions and outlook.

2 Numerical modeling of drop-freezing experiments

2.1 HUB-forward method to compute the fraction of frozen droplets f_ice(T) and cumulative freezing spectrum N_m(T) from a known underlying distribution P_u(T)

In the HUB-forward analysis we know or assume an underlying distribution P_u(T) of ice nucleation temperatures for the IN in the sample and generate from it an artificial IN dilution series similar to those obtained in experiments, from which we compute the cumulative freezing spectrum N_m(T) using Vali's equation (Eq. 1a). Using this approach, we investigate the relationship between N_m(T) and P_u(T) (Fig. 1) and the sensitivity of N_m(T) with respect to the number of droplets and dilutions. For generality, we represent P_u(T) as a linear combination of normalized continuous distributions P_i(T) that represent subpopulations of freezing temperatures:

$$\begin{array}{}\text{(2)}& P_{\mathrm{u}}\left(T\right)=c_{\mathrm{1}}{P}_{\mathrm{1}}\left(T\right)+c_{\mathrm{2}}{P}_{\mathrm{2}}\left(T\right)+\mathrm{\dots }+c_p{P}_p\left(T\right),\end{array}$$

where p is the total number of subpopulations, $P_{\mathrm{1}}\left(T\right),P_{\mathrm{2}}\left(T\right),\mathrm{\dots },P_p\left(T\right)$ are normalized distribution functions, and c₁, c₂, …, c_p are their weights such that $\sum _{i=\mathrm{1}}^p{c}_i=\mathrm{1}$. These subpopulations could correspond to different chemical, topographical, or structural motifs in the IN samples, although chemically distinct species could also produce overlapped freezing signatures, and a single species could display a broad freezing range. Our formalism does not require a mapping of subpopulations of freezing temperatures to physical IN sites. The units of P_u(T) are the same as for n_m(T), i.e., those of the cumulative spectrum divided by a unit of temperature, but are generally omitted in what follows. Throughout this work we assume that P_i(T) can be represented by Gaussian (i.e., normal) distributions:

$$\begin{array}{}\text{(3)}& P_i\left(T\right)=\left({\displaystyle \frac{\mathrm{1}}{s_i\sqrt{\mathrm{2}\mathit{\pi }}}}\right)e^{-\frac{\mathrm{1}}{\mathrm{2}}{\left[\frac{T-T_{\mathrm{mode},i}}{s_i}\right]}^{\mathrm{2}}},\end{array}$$

where each subpopulation P_i(T) is further characterized by its most likely temperature of freezing T_mode,i and spread of distribution of freezing temperatures s_i. We also provide in the HUB code the option for the user to use the log-normal distribution, which has a tail towards higher temperatures, or the left-tailed Gumbel distribution, which has a tail towards lower temperatures. In our model, we assume that the underlying distribution of ice-nucleating temperatures P_u(T) does not change with the concentration of INs. This last condition is violated when INs are involved in chemical, aggregation, or solubility equilibria that alter the proportionality between their concentration and the dilution factor of the sample, resulting in a lack of overlap of the pieces of the cumulative spectra N_m(T) obtained from different dilutions (Bogler and Borduas-Dedekind, 2020).

The number of INs in each droplet is given by the Poisson distribution:

$$\begin{array}{}\text{(4)}& p(n,\mathit{\lambda })={\displaystyle \frac{{\mathit{\lambda }}^n}{n\mathrm{!}}}{e}^{-\mathit{\lambda }},\end{array}$$

where n is the actual number of INs in each droplet and λ represents the average number of INs among all droplets of the corresponding dilution. Figure 2a shows the probability mass function (PMF) for λ=1, 5, and 10, computed according to Eq. (4) and sampling over N₀=10⁴ droplets using the “SciPy Stats” Python framework (Virtanen et al., 2020). As λ increases, the probability that any droplet nucleates homogeneously rapidly approaches zero (inset of Fig. 2a). When there is one IN on average per droplet $(\mathit{\lambda }=\mathrm{1}),\sim \mathrm{37}$ % of the droplets do not have any INs; i.e., they are “empty” droplets that would nucleate at the homogeneous nucleation temperature. We note that by performing dilutions until a sizeable fraction of droplets nucleate homogeneously, it is possible to calibrate the absolute concentration of ice nuclei in the original, undiluted sample.

To illustrate how the heterogeneous ice nucleation temperatures recorded in drop-freezing experiments depend on the number of INs in the droplets, we start from two examples with P_u(T) represented by one or two Gaussian subpopulations, shown with dashed black lines in Fig. 2b and c, respectively. We assign a temperature to each IN contained in droplets from a 10-fold dilution series of five solutions with λ=1, 10, 10², 10³, and 10⁴ average number of INs per droplet. If the droplet volume is constant, λ is proportional to the concentration of INs in the droplets. We sample N₀=10⁴ droplets for each concentration. This N₀ is much higher than the ∼100 droplets usually sampled in laboratory experiments; we address the effect of sampling in Sect. 3.1 below.

To sample independent random values for each IN, the number of random variates, which are drawn from P_u(T), is the total number of INs among N₀ droplets. Thereby, each droplet has a set of temperatures $T_j^{\mathit{\lambda }}=(T_{\mathrm{1}}^{\mathit{\lambda }},T_{\mathrm{2}}^{\mathit{\lambda }},\mathrm{\dots },T_k^{\mathit{\lambda }})$, where j is the droplet index and k is the IN index. Since we assume that freezing occurs at the characteristic temperature of the IN with the highest freezing temperature, the nucleation temperature for each droplet is defined as the maximum, i.e., the extreme upper value, of several independent freezing temperatures $T_{\mathrm{het},j}^{\mathit{\lambda }}=max(T_{\mathrm{1}}^{\mathit{\lambda }},T_{\mathrm{2}}^{\mathit{\lambda }},\mathrm{\dots },T_k^{\mathit{\lambda }})$. Figure 2b and c show the normalized distribution of $T_{\mathrm{het},j}^{\mathit{\lambda }}$ for different values of λ, namely $P_{\max }^{\mathit{\lambda }}\left(T\right)$. Therefore, P_u(T) represents the underlying probability of heterogeneous ice nucleation temperatures independent of the concentration of INs, while $P_{\max }^{\mathit{\lambda }}\left(T\right)$ represents the concentration-dependent distribution and has the same units as P_u(T) and the differential spectrum. According to the Fisher–Tippett–Gnedenko theorem, the distribution of extreme upper values of the Gaussian distribution is the right-skewed Gumbel distribution (Castillo et al., 2005; David and Nagaraja, 2004; Gumbel, 2012; De Haan and Ferreira, 2006), which has a fatter tail on the high-temperature side of its maximum. Indeed, the shift in $P_{\max }^{\mathit{\lambda }}\left(T\right)$ curves in Fig. 1b and c evinces that as the number of INs in the droplet increases, the probability of sampling the higher temperature tail of P_u(T) increases significantly. This skew is the reason why several dilutions are needed to sample the full population of ice nucleants.

Figure 2 (a) Probability mass function (PMF) of the Poisson distribution representing the number of INs per droplet. Colors represent different average numbers of INs per droplet: λ=1 (blue squares), λ=5 (purple triangles), and λ=10 (cyan circles). The inset shows the fraction of empty droplets as a function of λ. The connecting lines are solely guides for the eye. Panels (b) and (c) show the normalized underlying distributions P_u(T) of heterogeneous ice nucleation temperatures (dashed magenta line), composed of two subpopulations and one subpopulation, respectively. Colors represent the concentration-dependent normalized distribution $P_{\max }^{\mathit{\lambda }}\left(T\right)$ of heterogeneous ice nucleation temperatures: λ=1 (blue squares), λ=10 (cyan circles), λ=10² (green diamonds), λ=10³ (yellow ×), and λ=10⁴ (red triangles) INs per droplet. A bin width of 0.1 was used for P_u(T) and $P_{\max }^{\mathit{\lambda }}\left(T\right)$. All distributions were obtained using 10⁴ droplets. While the HUB-forward code explicitly accounts for $N_{\mathrm{F}}^{\mathit{\lambda }}$ and N₀, we note that their ratio can be approximated by $N_{\mathrm{F}}^{\mathit{\lambda }}/N_{\mathrm{0}}\approx (\mathrm{1}-e^{-\mathit{\lambda }})$ based on properties of the Poisson distribution.

HUB-forward computes the fraction of frozen droplets and cumulative spectra from a proposed underlying distribution of freezing temperatures, using extreme value statistics. The fraction of frozen droplets $f_{\mathrm{ice}}^{\mathit{\lambda }}\left(T\right)$ can be calculated as a function of the concentration-dependent distribution:

$$\begin{array}{}\text{(5)}& f_{\mathrm{ice}}^{\mathit{\lambda }}\left(T\right)=\underset{T_{\mathrm{m}}}{\overset{T}{\int }}{P}_{\max }^{\mathit{\lambda }}\left(T^{\prime }\right)\mathrm{d}{T}^{\prime }\times {\displaystyle \frac{N_{\mathrm{F}}^{\mathit{\lambda }}}{N_{\mathrm{0}}}},\end{array}$$

where the integration is from the ice melting temperature T_m to the temperature $T,N_{\mathrm{F}}^{\mathit{\lambda }}$ is the total number of droplets that freeze heterogeneously, and N₀ is the total number of droplets. We note that the approach taken in this work differs from that of previous studies that start from a microscopic model for the nucleation sites and nucleation theory to predict the fraction of frozen droplets using Monte Carlo simulations, as well as from previous modeling using the singular approximation, which do not account for the statistics of extreme sampling.

To use the HUB-forward code, the user must define the total number of droplets “ndroplets” that serves as the total number of each concentration and the number of subpopulations “nsubpop”. If “nsubpop” = 1, the user must provide the temperature of maximum likelihood T_mode,1 and the spread s₁. If “nsubpop” = 2, the user must provide T_mode,1, s₁, T_mode,2, s₂ and c₂. If “nsubpop” = 3, the user has to provide T_mode,1, s₁, T_mode,2, s₂, c₂, and T_mode,3, s₃, c₃. To generate the cumulative freezing spectrum N_m(T), the user needs to define the total number of concentrations “nconc”; the concentration of the parent suspension is defined in “density” and the droplet volume in “volumedrop”. The output is composed of different data plots and files: the normalized P_u(T) and $P_{\max }^{\mathit{\lambda }}\left(T\right)$, the artificially generated $f_{\mathrm{ice}}^{\mathit{\lambda }}\left(T\right)$, and N_m(T) built from the 10-fold dilution series.

Figure 3a and b show the fraction of ice computed using $P_{\max }^{\mathit{\lambda }}\left(T\right)$ of Fig. 2b and c, which correspond to P_u(T) with two subpopulations and one subpopulation, respectively. The intermediate plateau in Fig. 3b indicates that no droplets freeze at those temperatures. As discussed above, only 63 % of the droplets freeze heterogeneously for λ=1. We assume droplets of uniform volume V_drop=0.1 µL obtained through 10-fold dilution of a parent suspension with λ=10⁴ INs per droplet corresponding to a mass m=1 mg of IN in a volume V_wash=1 mL. We use Eq. (5) and the P_u(T) of Fig. 2b and c to generate $f_{\mathrm{ice}}^{\mathit{\lambda }}\left(T\right)$ (Fig. 3a and b), sampling either 100 or 10⁴ droplets per dilution. We combine the $f_{\mathrm{ice}}^{\mathit{\lambda }}\left(T\right)$ using Eq. (1a) to build the cumulative freezing spectra N_m(T) shown in Fig. 3c and d (sampling 10⁴ droplets per dilution) and Fig. 3e and f (sampling 10² droplets per dilution).

Figure 3 Panels (a) and (b) represent the fraction of ice $f_{\mathrm{ice}}^{\mathit{\lambda }}\left(T\right)$ computed using Eq. (5) and artificially generated data using 10⁴ droplets. Panels (c) and (d) are the corresponding cumulative freezing spectra N_m(T) computed using Vali's equation. Colors represent the different number of INs per droplet: λ=1 (blue squares), λ=10 (cyan circles), λ=10² (green diamonds), λ=10³ (yellow ×), and λ=10⁴ (red triangles). Panels (e) and (f) represent N_m(T)  obtained using 100 droplets. The dashed black lines in (c) and (d) indicate the temperatures corresponding to the location of the mode(s) in the underlying distribution. The dotted magenta lines are the knee points computed with the Python function “kneed”.

The ability of HUB-forward to generate the cumulative freezing spectrum N_m(T) from the underlying distribution P_u(T) allows for an analysis of the sensitivity of N_m(T) and P_u(T) to the number of droplets and dilutions, as seen in the comparison of N_m(T) generated from the same underlying distributions using 100 and 10⁴ droplets in Fig. 3. In Sect. 3.1 we show that the sampling with 100 droplets for only four dilutions of a system with two subpopulations of INs results in distortions of the distribution of freezing temperatures and the proportions of these populations in the differential spectrum.

The knee point in N_m(T) corresponds to the point of maximum curvature (Satopaa et al., 2011) and has been used to characterize the nucleation temperature of a particular subpopulation (Hartmann et al., 2022). Similar to Hartmann et al. (2022), we have identified in Fig. 3c and d the knee points (dotted magenta line) of the artificially generated N_m(T) by using a Python function named “kneed”. The Python function “kneed” using S=1, curve = “concave”, and direction = “decreasing”. The knee points T_knee are very close to the temperatures of maximum likelihood T_mode (dashed black lines) of the corresponding underlying distribution P_u(T), because under these conditions the differential freezing spectrum n_m(T) is a very good approximant for P_u(T). However, we find that the removal of the more dilute solutions eliminates the plateau in N_m(T) and results in poor estimation of the modes of P_u(T) from the knee of N_m(T).

2.2 HUB-backward method to recover the differential freezing spectrum n_m(T) from the cumulative freezing spectrum N_m(T) by a stochastic optimization procedure

The HUB-backward code implements a stochastic optimization procedure to extract the differential spectrum n_m(T) from a given cumulative spectrum N_m(T) or from an experimental f_ice(T) curve. The latter is useful when data are available for a single concentration. One possibility for obtaining n_m(T) from N_m(T) would be to follow the following steps: (i) propose a trial function $n_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$, (ii) use HUB-forward to predict the concentration-dependent distributions $P_{\max }^{\mathit{\lambda },\mathrm{trial}}\left(T\right)$ for various IN concentrations, (iii) use these in Eq. (5) to predict the freezing fractions $f_{\mathrm{ice}}^{\mathit{\lambda },\mathrm{trial}}\left(T\right)$, (iv) compute $N_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$ from the freezing fractions using Eq. (1a), (v) evaluate the difference between that trial and the target (experimental) value using

$$\begin{array}{}\text{(6)}& \mathit{\delta }\left(T\right)=\left|{\log }_{\mathrm{10}}\left[N_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)\right]-\left.{\log }_{\mathrm{10}}\left[N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)\right]\right|,\right.\end{array}$$

and then (vi) evolve the parameters that determine $n_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$ until the difference δ(T)  is minimized. However, the use of HUB-forward in steps (ii) and (iii) to generate and evaluate hundreds of droplets containing up to tens of millions of INs would require significant computations that render this optimization process inefficient.

The HUB-backward optimization procedure, sketched in Fig. 4, uses a shortcut for steps (ii) and (iii) above to directly predict $N_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$ from $n_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$ with fast convergence. The shortcut is based on the understanding that, in the asymptotic limit in which the sample is extremely dilute (i.e., λ⟶0), each droplet that nucleates heterogeneously contains a single IN. In such a case, sampling an infinitely large number of droplets with $P_{\max }^{\mathit{\lambda }⟶\mathrm{0}}\left(T\right)$ is equivalent to sampling each and every IN, i.e., $P_{\max }^{\mathit{\lambda }⟶\mathrm{0}}\left(T\right)=P_{\mathrm{u}}\left(T\right)$. In agreement with this ansatz, Fig. 2b and c show that the underlying distribution P_u(T) (dashed black line) and the concentration-dependent $P_{\max }^{\mathit{\lambda }=\mathrm{1}}\left(T\right)$ (blue squares) sampled with 10⁴ droplets per dilution are already very close, i.e., $P_{\mathrm{u}}\left(T\right)\approx P_{\max }^{\mathit{\lambda }=\mathrm{1}}\left(T\right)$.

Figure 4 Flowchart of the optimization procedure to obtain the differential freezing spectrum n_m(T) from the full cumulative freezing spectrum N_m(T) or fraction of frozen droplets f_ice(T).

With this insight and considering the intrinsic cumulative spectrum, $I_{\mathrm{u}}\left(T\right)={\int }_{T_{\mathrm{m}}}^T{P}_{\mathrm{u}}\left(T^{\prime }\right)\mathrm{d}{T}^{\prime }\times (\mathrm{1}-e^{-\mathrm{1}})$, we define the cumulative integral of the differential spectrum as

$$\begin{array}{}\text{(7)}& I_{\mathrm{u}}^{\mathrm{trial}}\left(T\right)=\underset{T_{\mathrm{m}}}{\overset{T}{\int }}{n}_{\mathrm{m}}^{\mathrm{trial}}\left(T^{\prime }\right)\mathrm{d}{T}^{\prime }\times \mathit{\beta },\end{array}$$

where the integration is from the ice melting temperature T_m to the temperature T and β is an adjustable scaling factor to be obtained from the optimization. Likewise, a similar estimate can be made for a single fraction of ice curve $f_{\mathrm{ice}}^{\mathrm{trial}}\left(T\right)=I_{\mathrm{u}}^{\mathrm{trial}}\left(T\right)$ using Eq. (7) and the mean squared error can be directly evaluated (Fig. 4). When the target is a cumulative freezing spectrum, HUB-forward uses $n_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$ to predict a trial cumulative freezing spectrum (Fig. 4),

$$\begin{array}{}\text{(8)}& N_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)=-\ln [\mathrm{1}-I_{\mathrm{u}}\left(T\right)]\times {\displaystyle \frac{\mathrm{1}}{X}},\end{array}$$

where $\mathrm{1}/X$ corresponds to the maximum of the cumulative spectrum in the target distribution, $\mathrm{1}/X=max\left[N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)\right]$. With Eq. (8) we obtain an $N_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$ that we compare with the target using Eq. (6) (Fig. 4). To do the comparison, HUB-backward uses a spline fit to interpolate the experimental $N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)$, in order to have equally spaced temperature points to compare with the estimates in $N_{\mathrm{m}}^{\mathrm{trial}}\left(T\right)$. We use the “interp1d” algorithm, which is available in the Python SciPy library (Virtanen et al., 2020), with a linear interpolation to construct new equally spaced data points within the range of the lowest and highest temperature values in the freezing spectrum. The cost function for the optimization is the mean squared error (MSE), computed from the difference δ(T) in Eq. (6):

$$\begin{array}{}\text{(9)}& \text{MSE}={\displaystyle \frac{\mathrm{1}}{t}}\sum {\mathit{\delta }}^{\mathrm{2}},\end{array}$$

where t represents the total number of equally spaced points in δ(T).

We use a stochastic global optimization technique based on a simulated annealing algorithm to find the set of parameters of $n_{\mathrm{m}}^{\mathrm{trial}}$ (Eqs. 2 and 3) and β (Eq. 7) that globally minimize the MSE. We use the simulated annealing (SA) algorithm “dual annealing” that is part of the SciPy minimize library (Virtanen et al., 2020) with its default arguments predefined, except for the parameters “maxfun” that sets the maximum number of evaluations of the objective function (we select “maxfun” = 1 000 000 in the examples below) and the seed for the generation of random numbers (a new random integer is automatically generated every time the HUB-backward code is run). We show below that the optimized differential spectra, $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$, are quite insensitive to the value of the seed.

The output of HUB-backward is an optimized differential spectrum $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ or an optimized fraction of ice $f_{\mathrm{ice}}^{\mathrm{optimized}}\left(T\right)$. To quantify how much this optimized prediction deviates from the known underlying distribution in the examples of Fig. 5, where P_u(T) is known, we define the mean relative error (MRE) for the set of parameters:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\text{MRE}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{3}p}}\sum _{i=\mathrm{1}}^p\left[\left|{\displaystyle \frac{T_{\mathrm{mode},i}^{\mathrm{optimized}}-T_{\mathrm{mode},i}^{\mathrm{target}}}{T_{\mathrm{mode},i}^{\mathrm{target}}}}\right|\right.\\ & \left.+\left|{\displaystyle \frac{s_i^{\mathrm{optimized}}-s_i^{\mathrm{target}}}{s_i^{\mathrm{target}}}}\right|+\left|{\displaystyle \frac{c_i^{\mathrm{optimized}}-c_i^{\mathrm{target}}}{c_i^{\mathrm{target}}}}\right|\right],\end{array}\end{array}$$

where p is the number of subpopulations.

We now turn our focus to how to select the input parameters required by HUB-backward to start the search for the underlying distribution, using the experimental $N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)$ or $f_{\mathrm{ice}}^{\mathrm{target}}\left(T\right)$ as a guide. The code requires the user to define the number of distinct Gaussian subpopulations P_i(T) that comprise the underlying distribution (Eq. 2) and to provide upper and lower bounds for the weights c_i, their modes T_mode,i, and spreads s_i of each of these populations. In general, we find that defining the minimum and maximum values for the weights to $c_i^{\max }=\mathrm{1}$ and $c_i^{\min }=\mathrm{0}$ (see constraint in Eq. 2), for the modes $T_{\mathrm{mode},i}^{\max }$ and $T_{\mathrm{mode},i}^{\min }$ to between the homogeneous nucleation temperature (about −30 ^∘C) and the melting temperature (0 ^∘C), and for the spreads to $s_i^{\max }=\mathrm{10}$ ^∘C and $s_i^{\min }=\mathrm{0.1}$ ^∘C works well. However, these bounds can be tuned in order to better fit the data (as we find for pollen in Sect. 3.2 below). If the existing experimental $N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)$ data are very noisy, they can be interpolated in HUB-backward using the method “interp1d” with “npoints” = 100 and then smoothed with a Savitzky–Golay filter by changing the parameters “window_length”, which is the length of the filter window, and “polyorder”, which is the order of the polynomial used to fit the samples (“filter” in Fig. 4). The default values are 3 and 1, respectively. HUB-backward generates a plot that compares the original and the interpolated target data.

To identify the minimum number of subpopulations needed to represent a given freezing spectrum, we consider that every time a population is accumulated in N_m(T) or f_ice(T), these functions display a sharp increase. We note that assuming a large number of subpopulations may challenge the interpretability of the optimized differential spectrum $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$.

We apply the HUB-backward procedure to the N_m(T) obtained in Fig. 3c and d by sampling four 10-fold dilutions with 100 droplets, i.e., only a total of 500 droplets. Figure 5 shows the comparison between the predicted (solid magenta lines) and the target (dashed black lines) N_m(T) (panels a and b) and n_m(T) (panels c and d). Table 1 shows the predicted parameters and the precision of the optimization procedure to recover the known underlying distribution P_u(T). The MRE between the underlying distribution P_u(T) and the optimized differential spectrum $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ is 2 % for the system with one subpopulation and 13 % for the one with two despite the low number of droplets used to sample the cumulative freezing spectra in the computer-generated freezing experiments.

Table 1 Mean relative error (MRE), mean squared error (MSE), and parameters of the optimized differential freezing spectra $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ obtained using the HUB-backward code. The values shown here were calculated based on the average of n=3 independent runs. The error bars, shown in parentheses, were calculated by dividing the standard deviation of the values in these runs by 3${}^{\mathrm{1}/\mathrm{2}}$.

Figure 5 Panels (a) and (b) show the comparison between $N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)$ (black circles) and the $N_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ computed with the optimized solution $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ (solid red line). Panels (c) and (d) show the known underlying distributions P_u(T) (dashed black line) and the optimized underlying distributions $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ (solid red line) based on three independent runs. The parameters of the predicted underlying distribution $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ are summarized in Table 1.

We conclude that the HUB-backward code gives a good estimate of the mode, spread, and weights of the populations of INs in a sample and it can be applied in a situation where P_u(T) is unknown. In Sect. 3.1 we discuss how the accuracy is of the underlying distribution recovered with HUB-backward impacted by various schemes of the sampling of the number of droplets and dilutions to construct N_m(T). In Sect. 3.2, we apply the HUB-backward procedure to obtain $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ from actual N_m(T) of experiments with various soluble biological INs. In Sect. 3.3, we apply the HUB-backward procedure to obtain $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ from f_ice(T) of experiments of insoluble crystal INs.

3 Using the HUB code to optimize and analyze drop-freezing experiments

3.1 Effect of the number of droplets and dilutions on the cumulative freezing spectrum N_m(T)

Figure 3d–f show N_m(T) generated with HUB-forward using five dilutions from λ=10⁴ to 1 of a solution with P_u(T) containing two populations in a ratio of 9 to 1. The N_m(T) are different when the number of droplets per dilution is 100 (Fig. 3f) or 10⁴ (Fig. 3d). As shown in the previous section, the freezing spectrum obtained with 100 droplets and five dilutions has enough sampling to recover this P_u(T) with good accuracy (Fig. 5c and d). We test different number of droplets and concentrations, defined by the average number of INs per droplet λ, to test the sensitivity of n_m(T) to the number of droplets and dilutions when the underlying distribution P_u(T) is known. We use HUB-forward to build N_m(T) based on a combination of different numbers of droplets and concentrations, similar to the case shown in Fig. 3f. Then, we use HUB-backward to obtain $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$, compare it to P_u(T), and test the accuracy of each prediction through its mean relative error (MRE) defined in Eq. (10).

The left panels of Fig. 6 show N_m(T) generated with HUB-forward based on a combination of different concentrations using 100 droplets for each dilution. The magenta lines are based on the data provided by the HUB-backward code. The right panels of Fig. 6 compare $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ in magenta and the known underlying distribution P_u(T) in black. In this example, n_m(T) is very close to P_u(T) if both subpopulations are sampled enough. However, if the most dilute solution with λ=1 is not included in N_m(T) (second panel), the estimate of the underlying distribution is very poor. Thus, to improve the sampling of the lower tail of P_u(T), we recommend ending the dilution series always in the immediacy of λ=1, which can be gleaned from the temperature range for which N_m(T) becomes flat and a sizeable fraction of droplets of the more diluted sample nucleates homogeneously (inset of Fig. 2a). The one-to-one correspondence between the fraction of droplets nucleated homogeneously and the average number of particles in the droplet in the highly diluted limit (inset of Fig. 2a) demonstrates that reaching this limit allows for an absolute calibration of the number of INs in the initial sample. Moreover, sampling to concentrations down to about one nucleant per droplet is essential to recover a proper weight of the poorly nucleating IN populations.

Figure 6 Panels (a, c, e, g) represent the cumulative freezing spectra N_m(T) sampled from the same underlying distribution P_u(T). Colors represent the different numbers of INs per droplet: λ=1 (blue squares), λ=10 (cyan circles), λ=10² (green diamonds), λ=10³ (yellow ×), and λ=10⁴ (red triangles). The sampling was done using 100 droplets for each concentration. Panels (b, d, f, h) represent the differential freezing spectra n_m(T) compared to the known underlying distribution P_u(T), shown by the magenta and dashed black lines, respectively. Panels (a–d) were computed with a different number of dilutions. The mean relative error (MRE) was computed using Eq. (10). The parameters of n_m(T) and P_u(T) are shown in Table S1.

The relative weights of class A and C populations in Pseudomonas syringae is approximately 1 to 1000 (Sect. 3.2), while the ratio is 9 to 1 in the two-population system example of Fig. 6. To understand the impact of highly imbalanced populations on the sampling of the cumulative spectrum and recovery of the underlying distribution, we show in Fig. 7 the analysis of an example where the subpopulation of highly efficient INs is 3 orders of magnitude less likely to occur than the subpopulation at lower temperatures, mimicking the one of P. syringae. Our analysis confirms that it is important to end the dilution series in the immediacy of λ=1 to fully represent the contribution of the poorer INs (Fig. 7b–f). Furthermore, we find that it is important to sample a concentration high enough to account for the rare INs that nucleate at the warmest temperatures (Fig. 7d–h).

If only 25 droplets per dilution, instead of 100, are used to construct the cumulative spectrum, the impact of insufficient sampling at the higher concentrations is more pronounced: compare Fig. 8c and Fig. 7d obtained with the same underlying distribution P_u(T) with 1000 to 1 subpopulation ratios and number of dilutions.

We conclude that an increase in the accuracy in the account of the subpopulations requires a higher number of dilutions and the checking of the predictions with the addition of each successive concentration to ensure convergence of $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$. Measuring fewer droplets or fewer dilutions leads to poor statistics and results in incompleteness or the misrepresentation of the underlying distribution in samples with multiple subpopulations. In principle, increasing the number of droplets of the most concentrated solutions or adding more 10-fold concentrated ones until there are no changes in the cumulative spectrum is recommended to ensure complete sampling. When that limiting scenario is not attainable, the use of HUB-forward to produce synthetic data from a proposed underlying distribution, followed by the recovery of the differential spectrum from these data sets, allows for an estimation of the errors that may be incurred for putative, proposed underlying distributions with the sampling scheme available in the laboratory.

Figure 7 Panels (a–d) represent the cumulative freezing spectra N_m(T) sampled from the same underlying distribution P_u(T). Colors represent the different number of INs per droplet: λ=1 (blue squares), λ=10 (cyan circles), λ=10² (green diamonds), λ=10³ (yellow ×), and λ=10⁴ (red triangles). The sampling was done using 100 droplets for each concentration. Panels (e–h) represent the differential freezing spectra n_m(T) compared to the known underlying distribution P_u(T), shown by the magenta and dashed black lines, respectively. The mean relative error (MRE) was computed using Eq. (10) and the parameters of n_m(T) and P_u(T) are shown in Table S2.

Figure 8 Panels (a–c) represent the cumulative freezing spectra N_m(T) sampled from the same underlying distribution P_u(T). Colors represent the different number of INs per droplet: λ=1 (blue squares), λ=10 (cyan circles), λ=10² (green diamonds), λ=10³ (yellow ×), and λ=10⁴ (red triangles). The sampling was done using 25 droplets for each concentration. Panels (d–f) represent the differential freezing spectra n_m(T) compared to the known underlying distribution P_u(T), shown by the magenta and dashed black lines, respectively. The mean relative error (MRE) was computed using Eq. (10). The parameters of n_m(T) and P_u(T) are shown in Table S3.

3.2 Obtaining the differential freezing spectrum n_m(T) from the experimental cumulative freezing spectrum N_m(T) of biological INs using the HUB-backward code

In this section we use the HUB-backward code to obtain the differential freezing spectrum n_m(T) from the cumulative freezing spectra N_m(T) of the fungi Fusarium acuminatum strain 3–68 (Kunert et al., 2019), the bacterium P. syringae (Schwidetzky et al., 2021), and birch pollen (Dreischmeier, 2019). We select these systems because they are important biological INs and show increasing complexity in terms of the apparent number of underlying distributions that define their freezing spectra.

The experimental N_m(T) obtained for F. acuminatum (black squares in Fig. 9a) was obtained by sampling six 10-fold dilutions, each with 96 droplets (Kunert et al., 2019). Figure 9a shows the cumulative spectra optimized assuming one (green curve) and two (cyan curve) subpopulations; Fig. 9b shows the corresponding optimized differential freezing spectra. The $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ with a single subpopulation that peaks at −5.9 ^∘C is unable to represent the cumulative density of the most potent nuclei and misses the inflection at around −5.9 ^∘C in the experimental data, resulting in a mean squared error MSE = 0.05. The $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ with two subpopulations has a lower MSE = 0.003 and a better fit that suggests a population that peaks at −7.3 ^∘C and another at −5.5 ^∘C, in comparable amounts (Table 2). Most notably, the two subpopulations do not overlap in the differential freezing spectrum, supporting that they may indeed correspond to different physical entities. The improvement in the fit becomes apparent in the inset of Fig. 9a, which shows N_m(T) on a linear scale. The significant slope of N_m(T) even at the lowest temperatures indicates that the sampling of more diluted solutions is needed to capture the contribution of the less active INs. An attempt to represent F. acuminatum nucleation data with three different subpopulations resulted in two of them being almost identical. We conclude that adding a third subpopulation is unnecessary to reproduce the experimental cumulative freezing spectrum of F. acuminatum. We refer the reader to Schwidetzky et al. (2023) for an interpretation of the size of the ice-nucleating surface of F. acuminatum based on its differential spectrum and nucleation theory.

Figure 9 Cumulative freezing spectra N_m(T) obtained from drop-freezing experiments for (a) F. acuminatum strain 3–68 (Kunert et al., 2019), (b) P. Syringae (Schwidetzky et al., 2021), and (c) birch pollen (Dreischmeier, 2019) (black circles). The solid green, long dashed cyan, and short dashed red lines represent $N_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ computed with the optimized differential freezing spectra $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ obtained with the HUB-backward code considering one, two, and three subpopulations, respectively. Panels (b), (d), and (e) show $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$. The gray circles are experimental data points in the measurement of the birch pollen ice nucleation spectrum that were not considered in the optimization procedure. Insets in (a) and (c) show N_m(T) in normal scale.

Table 2 Mean squared error (MSE) and parameters of the differential freezing spectra n_m(T) obtained using the HUB-backward code and experimental data as input. The values shown here were calculated based on the average of n=3 independent runs. The error bars, shown in parentheses, were calculated by dividing the standard deviation of the values in these runs by 3${}^{\mathrm{1}/\mathrm{2}}$.

Next, we apply the HUB-backward code to analyze the experimental freezing spectrum of Snomax^®, i.e., inactivated P. syringae. The cumulative spectrum suggests the presence of two distinct subpopulations, usually called class A (at warmer temperatures) and class C (at colder ones). We first assume the differential freezing spectrum n_m(T) of P. syringae is a combination of two Gaussian populations. The parameters of the optimized differential spectrum with two subpopulations are listed in Table 2, and the curve is shown in Fig. 9d with a cyan line. We use a logarithmic scale to represent this $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ because the population corresponding to class A accounts for less than 0.1 % of the total (Table 2). While the fit with two subpopulations results in a good overall account of the target data, we note that there is some difference in the region between classes A and C (Fig. 9c). The fitting for P. syringae achieves an excellent agreement between optimized and target cumulative spectra (Fig. 9c), through the prediction of an additional peak located between classes A and C (the elusive class B), with a population comparable to class A (Table 2 and red curve in Fig. 9d). However, more measurements and analyses are needed to establish whether this “class B” peak at −5.2 ^∘C is reproducible and truly distinct from the one of class A at −3.7 ^∘C to warrant a physical interpretation. Overall, both the analyses with two and three subpopulations agree with previous ones (Govindarajan and Lindow, 1988; Warren, 1987) that concluded that over 99 % of the INs active in P. syringae bacteria in Snomax^® belongs to class C. The analysis presented here for fungal and bacterial INs illustrates how HUB-backward can be used to reveal and characterize the underlying number of IN subpopulations of complex biological samples.

To further test the methodology, we model the cumulative freezing spectrum of birch pollen. Given that the original N_m(T) data for pollen in Fig. 3.1 of Dreischmeier (2019) consist of multiple independent curves, we took one of the many presented in this graph as the target $N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)$ (black curve in Fig. 9e) and present some of the additional data – not used in the optimization – with gray circles in Fig. 9e. Section S4 in the Supplement shows that the differential spectrum optimized from the whole data set and its sparse sampling are almost identical because HUB-forward interpolates and smooths the input data to produce an equispaced data set. The $N_{\mathrm{m}}^{\mathrm{target}}\left(T\right)$ seems to contain three quite separated subpopulations, which is confirmed by the accuracy of the optimized cumulative spectrum in Fig. 9e. The parameters of the optimized differential freezing spectrum $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ and the MSE are shown in Table 2. Our analysis indicates that the two subpopulations that nucleate ice above −16 ^∘C constitute less than 0.01 % of the active nucleating sites in pollen (Fig. 9e), consistent with drop-freezing assays that only measured solutions with low concentrations of birch pollen and did not observe freezing at higher temperatures (Augustin et al., 2013; Pummer et al., 2012; Felgitsch et al., 2018), while the more extensive data of Dreischmeier (2019) reveal two more active subpopulations of INs.

To further illustrate the use of HUB-backward, Fig. 10 shows the effect of pH in the modes, spread, and weights of the subpopulations that contribute to the nucleation spectrum of P. syringae (Snomax^®), using data from Lukas et al. (2020). Freezing in the temperature range of class A drops about 3 orders of magnitude when the pH is lowered from 6.2 to 4.4 (Fig. 10b). However, we note that the cumulative number of INs is preserved in the experimental cumulative freezing spectrum (Lukas et al., 2020), indicating that the change in pH did not impact the number of nucleants. Figure 10c and d demonstrate that the distributions associated with both subpopulations shift to lower temperatures when the pH decreases, and the range of freezing temperatures in class A becomes broader. An attempt to fit the cumulative spectra of Snomax at different pH values with the same subpopulations, allowing only for adjustment of their weights, resulted in a poor fit to the experimental N_m(T), supporting the conclusions of Lukas et al. (2020) of a central role of electrostatic interactions in the assembly of the bacterial ice-nucleating proteins and their ability to bind to ice. This analysis exemplifies how HUB-backward can be applied to quantify the dependence of IN on environmental variables.

Figure 10 Effect of changing the pH on the subpopulations of P. Syringae (Lukas et al., 2020). (a) Differential freezing spectra n_m(T) obtained using the HUB-backward code. Colors represent the different pH values: 6.5 (long dashed black line), 5.6 (short dotted blue line), and 4.4 (solid magenta line). (b) Ratio between the weights, (c) the modes, and (d) the spreads of each subpopulation as a function of pH. The fitting of N_m(T) and the parameters of n_m(T) are shown in Fig. S1 and Table S4.

3.3 Obtaining the differential freezing spectrum n_m(T) from the experimental fraction of ice f_ice(T) of insoluble ice nucleators using the HUB-backward code

Section 3.1 and 3.2 discuss how to obtain the differential spectrum from a target cumulative one. However, there are many cases where the results are presented as a fraction of frozen droplets as a function of temperature, f_ice(T). In these cases, the HUB-backward code can be used to obtain the optimized differential freezing spectrum $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ directly from $f_{\mathrm{ice}}^{\mathrm{target}}\left(T\right)$. Section S5 illustrates this approach for the analysis of droplet freezing data for a sample of lignin (Bogler and Borduas-Dedekind, 2020) in which the INs participate in aggregation equilibria. Here, we exemplify the optimization of the differential spectrum of cholesterol from experimental freezing data obtained at two cooling rates with droplets sampled from a single dilution.

In the analysis of drop-freezing experiments, it is assumed that each IN has a singular freezing temperature, independent of the cooling rate. However, ice nucleation is a stochastic process, and the underlying distribution of freezing temperatures P_u(T) strictly depends on both temperature and cooling rate, as slower rates give more time for the system to cross the nucleation barrier at warmer temperatures.

The triangles and squares in Fig. 11a display the experimental f_ice(T) obtained by sampling the freezing of hundreds of 120 µL droplets pipetted from a suspension of cholesterol monohydrate crystals in contact with Teflon cooled at 0.18 K min⁻¹ (triangles) and 0.06 K min⁻¹ (squares) (Zhang and Maeda, 2022). Our analysis of the freezing data of cholesterol monohydrate shows that even a 3-fold change in the cooling rate can have a significant effect on the differential spectrum (Fig. 11b).

Figure 11 Use of the HUB-backward code to estimate the optimized differential freezing spectra $n_{\mathrm{m}}^{\mathrm{optimized}}\left(T\right)$ based on the fraction of frozen droplets $f_{\mathrm{ice}}^{\mathrm{target}}\left(T\right)$ of cholesterol (Zhang and Maeda, 2022) at different cooling rates. Black circles and squares are experimental data, and dashed cyan and solid red lines are the optimized differential spectra given by the HUB-backward code. The parameters of n_m(T) are shown in Table S5.

As expected, the modes of the three populations move towards warmer temperatures upon decreasing the cooling rate. We note, however, that the shift in the peaks is not uniform; the middle one seems to be more sensitive to the cooling rate. Different sensitivity of the freezing rate of subpopulations has also been reported in simulations of nucleation data of minerals using the stochastic and modified singular frameworks (Herbert et al., 2014; Murray et al., 2011). The modified singular model proposes an empirical correction to the relation between f_ice(T) and N_m(T) to account for the effect of the cooling rate on the shift in these quantities (Vali, 1994). That analysis could be extended to the analysis of the subpopulations of INs obtained with HUB-backward. Moreover, it would be interesting in future studies to use the rate dependence of the mode of the subpopulations to extract the steepness of the nucleation barrier with temperature using nucleation theory (Budke and Koop, 2015) and to investigate the relationship between the cooling rate dependence of the differential spectrum obtained in the singular approximation with the interpretation of the same data modeled with the stochastic framework, such as in Wright et al. (2013) and Herbert et al. (2014).

4 Conclusions

In this study, we present the HUB method and associated Python codes that model (HUB-forward code) and interpret (HUB-backward code) the results of droplet freezing experiments under the assumptions that each ice-nucleating site in the sample has a characteristic nucleation temperature that is time-independent. The use of the singular approximation is the same as that used by Vali (1971, 2014, 2019) in his derivation of the ice nucleation spectra from data of fraction of frozen droplets. Different to previous implementations of the singular model, HUB accounts for the distribution of the number of INs in droplets at a given concentration and uses extreme value statistics to represent the effect of dilutions in the frozen fraction and freezing spectra. Our method and codes allow users to obtain an analytical differential freezing spectrum n_m(T) from the experimental distribution of freezing temperatures, and vice versa. The differential freezing spectrum n_m(T) is an approximant to the underlying distribution of ice-nucleating temperatures P_u(T), which provides a hub to connect the experimental freezing temperatures with interpretative physical analyses using kinetic models or nucleation theory that can be used to elucidate the mechanisms of nucleation and origins of these distributions.

HUB-forward predicts the cumulative ice nucleation spectrum N_m(T) and fractions of frozen droplets f_ice(T) from a known (or assumed) underlying distribution P_u(T) of nucleation temperatures for the INs in the sample. The HUB-forward code can be used to investigate the effect of the number of droplets and dilutions on the temperature range of the cumulative freezing spectrum N_m(T). Our analysis shows that the differential freezing spectrum n_m(T) is identical to the underlying distribution of heterogeneous ice nucleation temperatures P_u(T) only when sampling is complete. Measuring fewer droplets or fewer dilutions can result in a biased representation of the differential and cumulative spectra. HUB-forward predicts f_ice(T) and N_m(T) from a proposed distribution of IN temperatures, allowing its users to test hypotheses regarding the role of subpopulations of nuclei in the freezing spectra and providing a guide for a more efficient collection of freezing data.

HUB-backward uses a non-linear optimization method to find the differential freezing spectrum n_m(T) that best represents the experimental target cumulative freezing spectrum N_m(T) or fraction of frozen droplets f_ice(T) in the experiments. The analytical form of the differential freezing spectrum n_m(T) obtained from HUB-backward offers an interpretable physical basis. The interpretability of the results in terms of subpopulations provides an advantage over polynomial fitting and differentiation of N_m(T). Indeed, we show that the HUB-backward code can be used to reveal and characterize the underlying number of IN subpopulations of complex biological samples (Snomax^®, fungi Fusarium acuminatum, and birch pollen) and quantify the dependence of their subpopulations on environmental variables. Interestingly, our analysis evinces subpopulations that are not obvious to the eye and have not previously been identified in these samples. The robustness of the signals that correspond to these populations and their physical nature require further investigation.

We illustrate the use of HUB-backward to obtain the differential freezing spectrum n_m(T) from the fraction of frozen droplets f_ice(T) collected at a single concentration. We apply that analysis to demonstrate that n_m(T) depends on the cooling rate. The shift in the peaks of the subpopulations to higher temperatures upon decreasing the cooling rate is not unexpected, as longer waiting times allow for the surmounting of the same nucleation barrier at warmer temperatures. By providing the temperature dependence of the mode, spread, and weight of the subpopulation peaks, HUB-backward can be combined with nucleation theory and other theoretical analyses to extract the steepness, and maybe even the distribution, of nucleation barriers that control the freezing process.