Multidimensional infrared diffusion-ordered spectroscopy in depletion mode distinguishes protein amyloids and monomers

The advent of two-dimensional infrared (2D-IR) spectroscopy, now precisely 25 years ago,¹ has improved our understanding of a wide range of molecular systems.^2–10 This holds, in particular, for our understanding of the formation and structure of amyloids.^11–19 Amyloid protein aggregates are widely studied because of their relevance in neurodegenerative diseases.²⁰ Parkinson’s disease, for instance, is caused by the aberrant aggregation of the protein α-synuclein. Such unwanted assembly leads to the formation of amyloid fibrils, which consist of repeats of β-sheet structures that can reach up to hundreds of nanometers in length.²¹ Amyloid formation is a multi-scale process: first, monomers aggregate to form oligomers, which then further assemble into fibrils. There is evidence that the oligomers, and not the mature fibrils, are the cytotoxic species in human diseases, increasing the urgency of investigating the formation and structure of amyloid aggregates at all relevant length scales.^22–24 

Infrared (IR) spectroscopy is very well suited for this purpose: It is a label-free, non-invasive method, and molecular vibrations, especially the amide I mode, are sensitive probes of protein structure.²⁵ This is because the couplings between the amide I vibrations in a protein depend on the spatial arrangement of the amide groups, giving rise to characteristic delocalized normal modes for different secondary structures.^25–31 In particular, β-sheets are characterized by two modes at 1620–1630 and 1680–1700 cm⁻¹,^{25,30,32–34} and the appearance of these bands in the FTIR spectrum, and their cross-peak pattern in the 2DIR spectrum, is a clear indication of the presence of amyloid oligomers³⁵ or fibrils.^{11–13,15–19,36} The application of conventional and two-dimensional IR spectroscopy has led to many new insights into amyloid structure and growth. However, it is not so easy to study mixtures of protein monomers, oligomers, and aggregates with IR-based methods since, in general, the amide I spectra of these species overlap strongly. As a consequence, kinetic studies of amyloid formation, in which the sample is inevitably a mixture of monomers, oligomers, and aggregates, often have to rely on curve-fitting strategies to disentangle the (2D)IR spectra.^18,37 Here, we show how this spectral-overlap problem can be solved using Infrared Diffusion-Ordered Spectroscopy (IR-DOSY)—a method that we recently developed³⁸ based on a similar concept in nuclear magnetic resonance (NMR) spectroscopy.^39–48 

IR-DOSY combines IR spectroscopy with microfluidic technology. The combination of these two methods has been already shown to be extremely effective to chemically image microfluidic flows in a label-free and noninvasive manner.^49–51 In situ chemical imaging is particularly useful in flow chemistry, in enabling the analysis of the fluid composition⁵⁰ or to monitor the evolution of reactions⁵¹ in real time within a microchannel. Recently, microfluidic technology was also combined with 2D-IR spectroscopy, and by exploiting the sensitivity of 2D-IR to molecular solvation and flow control of microfluidic devices, Tracy et al. were able to study the effect of solvent mixing on a test molecule.⁵² Here, we merge (2D)IR spectroscopy with microfluidic technology to achieve combined structure and size sensitivity. The information obtained is similar to that obtained by combining liquid chromatography with infrared spectroscopy, which involves somewhat more complicated setups.^53,54 Currently, the most commonly used methods for separating complex protein mixtures include flow-induced dispersion analysis (in which the diffusion constants of the components in a mixture are determined by observing the Taylor dispersion⁵⁵ in a hydrodynamic flow system⁵⁶) and tangential flow methods (in which the size sensitivity involves the use of membranes^57,58). We think that IR-DOSY could form a useful complement to these existing, more sophisticated methods.

In IR-DOSY, we use a relatively simple method to simultaneously measure the IR absorption and the diffusion coefficient of the molecules (or particles) in a sample, resulting in multidimensional spectra in which the diffusion coefficient is along one of the axes, and IR frequency (or frequencies, in the case of 2DIR) along the other. The idea is inspired by NMR-DOSY, but we achieve diffusion-coefficient selectivity in a different manner: instead of using a pulsed field gradient, as is done in NMR-DOSY, we create a spatially inhomogeneous distribution of the solute molecules inside an infrared sample cell and track the diffusion of the solute molecules in a IR-frequency resolved manner. Since the diffusion constant of a molecule is determined by its size through the Stokes–Einstein relation,⁵⁹ in an IR-DOSY spectrum of a mixture, the spectra of compounds with different sizes are cleanly separated into distinct sets of peaks at different positions along the diffusion-coefficient axis.³⁸ Since the protein monomers, oligomers, and amyloids have extremely different sizes, IR-DOSY seems ideally suited to unravel the (2D)IR spectra of their mixtures. The only practical problem is that the diffusion of amyloids is too slow to be experimentally observable. Here, we will show that with a small adaptation, IR-DOSY can still be used to separate the (2D)IR spectra of protein aggregates and monomers in a mixed sample. We will first test the idea on a mixed solution of two small molecules and then apply it to a mixture of bovine-serum-albumin amyloids and monomers.

Dialanine and acetone were dissolved in D₂O, to a final concentration of 40 mg/ml and 4% v/v respectively. Bovine serum albumin (BSA) aggregation was optimized from the procedure described by Militello et al.⁶⁰ A phosphate-buffered saline (PBS) solution of 0.1M was made by dissolving one tablet of PBS in D₂O (10.0 ml). The pD was measured to be 7.4 using a pH-meter (Thermo Scientific Orion 2-Star). The pH-values were corrected by +0.4 to obtain the pD-values. BSA was dissolved in PBS (0.1M) in an Eppendorf tube to various concentrations. The solution was shaken shortly with a Vortex mixer (Stuart) to ensure solvation of the protein. Next, the solution was centrifuged (Eppendorf Centrifuge 5418) at 2000 rpm for two minutes, and filtered through a 0.45 μm filter. Before filtering, the filter was wet with a bit of D₂O. The Eppendorf tube containing the protein solution was placed in a small glass vial in such a way that it did not touch the bottom of the vial. The vial was placed in the oven and incubated at different temperatures for varying times. Before and after heating, an IR-spectrum was measured, to check if fibrils were formed. The size of the fibrils was measured using dynamic light scattering.

For IR-measurements, a Bruker Vertex 70 was used. All measurements were performed at room temperature. Per measurement, 32 scans were made, with a spectral resolution of 2 cm⁻¹. The frequency range was selected from 7000 to 400 cm⁻¹. Infrared spectra were recorded every 60 seconds.

A detailed description of the setup used to measure the 2DIR spectra can be found in Ref. 61. Briefly, pulses of wavelength 800 nm and with a 40 fs duration are generated using a Ti:sapphire oscillator and further amplified using a Ti:sapphire regenerative amplifier to obtain 800 nm pulses at 1 kHz repetition rate. These pulses are converted in an optical parametric amplifier to obtain mid-IR pulses (∼20 μJ, ∼6100 nm) that have a spectral full width at half maximum (FWHM) of 150 cm⁻¹. The beam is split into a probe and a reference beam (each 5%), and a pump beam (90%) that is aligned through a Fabry–Pérot interferometer. The pump and probe beams are overlapped in the sample in an ∼250 μm focus. The transmitted spectra of the probe (T) and reference (T₀) beams with pump on and off are recorded after dispersion by an Oriel MS260i spectrograph (Newport, Irvine, CA) onto a 2 ×32 pixel mercury cadmium telluride (MCT) array. The probe spectrum is normalized to the reference spectrum to compensate for pulse-to-pulse energy fluctuations. The 2DIR signal is obtained by subtracting the probe absorptions in the presence and absence of the pump pulse.

For both IR-DOSY and 2DIR-DOSY experiments, an IR-DOSY cell with a spacer of 50 μm and a channel width of 4 mm was used, with a slit of 400 μm for the FTIR experiments. The IR-DOSY cell was flushed with solvent (about 0.25 ml) and, subsequently, with the sample (about 0.25 ml). A pump (Harvard PHD 2000 Infusion/Withdraw) was used to create a steady flow of both solutions. The flow rate was set to 25 μl/min.

To analyze the time-dependent spectral data and convert it into a DOSY spectrum, we use numerical solutions to the diffusion equation, which we described previously in Ref. 38, which we refer the reader to for the details. Briefly, we assume that N species are present in the solution, with IR (or equivalently the second-derivative) spectra A_i(ν) (with i the species number and ν the IR frequency). In depletion mode, the total absorption at IR frequency ν, position y = L/2, where L is the width of the sample-cell channel, and time t are given by

where D_i is the diffusion coefficient of species i, and C(L/2, D_it/L²) is defined as

with τ = Dt/L².

For τ ≪ 1, the sum converges slowly, and then it is more efficient to use the image-charge method,⁶² which gives

where $erf(x)=(2/π)∫0xe−t2dt$⁠. To efficiently calculate C(L/2, τ), we use Eq. (2) for τ > 0.2 and Eq. (3) for τ < 0.2. Using the first two terms of Eq. (2) and the terms −3 ≤ n ≤ 3 of Eq. (3), we obtain a precision at x = L/2 of $<10−8$ at all times (see the supplementary material).

By least-squares fitting Eq. (1), we obtain diffusion constants D_i and their associated spectra A_i(ν). The 2D-DOSY spectrum S(ν, D) is obtained by multiplying the spectral amplitude A_i(ν) with the appropriate probability distribution for D_i:⁴⁴ 

where N is the number of species (N = 2 in the experiments reported here), and σ_i is the uncertainty in the diffusion coefficient obtained from the least-squares fit.

Figure 1 shows schematically how we simultaneously determine the diffusion coefficients and IR spectra of the molecular species present in a mixed solution. The practical details of the setup and the data-analysis procedure can be found in Ref. 38. In a thin space (50–100 μm) between two IR-transparent windows, we inject the sample solution and its solvent with the same flow rates so that the contact area between the liquids is a line in the center of the cell (the flow is laminar³⁸). After we stop the flow, the solute molecules (blue and red dots in the figure) start to diffuse into the pure solvent region at a rate that depends on their diffusion coefficient (and, hence, on their size). Previously, we tracked the arrival of the diffusing molecules in the initially empty half of the sample volume,³⁸ but, here, we track the depletion of the IR signal in the mixture-filled half. This has the important advantage that the infrared spectrum of the species with very low diffusion constants is also detected (although their diffusion coefficient cannot be determined). Since the IR beam in a standard FTIR spectrometer typically has a diameter of several millimeters, we ensure that we measure the IR spectrum of this specific region of the sample-cell channel by means of an adjustable optical slit [green rectangle in Fig. 1(a)]. As time progress, the IR peaks of the solute molecules decrease as they diffuse into the empty half of the sample volume, at different rates depending on their size (in the figure: red first, blue later). The time dependence of the amplitude of each species is governed by the diffusion equation [see Eq. (1) above for the explicit solution], and from a global analysis of the time- and frequency-dependent data, we obtain an IR-DOSY spectrum in depletion mode [Fig. 1(b)]. Schematics of the cell and of the experimental implementation are shown in Figs. 1(c) and 1(d).

We first test IR-DOSY in depletion mode on a mixture of acetone and dialanine, which we previously investigated with IR-DOSY in “arrival” mode, monitoring the arrival of the diffusing species in the initially solvent-filled part of the sample cell.³⁸ Figure 2(a) shows the time-dependent IR-absorption spectrum recorded in the solution-filled part of the sample cell, at different times after the injection. The IR bands of acetone disappear faster than those of dialanine. In Fig. 2(b), we plot the absorption vs time at the frequency of the main absorption peaks of acetone and dialanine. The acetone signal decreases much faster, and the dialanine signal shows a short initial lag phase that is well reproduced by the diffusion equation (curves in the figure). By least-squares fitting the time- and frequency-dependent data with the solution of the diffusion equation (see the supplementary material for the details), we obtain the two diffusion coefficients and their associated species spectra. The diffusion coefficients obtained using IR-DOSY in depletion mode are similar to the ones found previously.³⁸ Figures 2(c) and 2(d) show the resulting IR-DOSY spectrum. The frequency region between 1350 and 1550 cm⁻¹ is generally highly congested because of the overlap of side-chain modes, such as the CH-bending modes. In the IR-DOSY spectrum in depletion mode, we observe that the spectral bands separate into two rows centered at D ∼ 3 · 10⁻⁶ and D ∼ 5.5 · 10⁻⁶ cm²/s. In the top row, the spectral bands are centered at 1695, 1450, and 1360 cm⁻¹, which correspond to the absorption frequencies of acetone (see the supplementary material for the spectrum of a solution of only acetone). In the bottom row, we observe spectral bands at 1665, 1600, 1490, 1400, and 1360 cm⁻¹, which correspond to the absorption frequencies of dialanine (see the supplementary material for the spectrum of dialanine solution). Just as for IR-DOSY in “arrival” mode, the IR-DOSY spectrum in depletion mode neatly resolves the congested FTIR spectrum [top panel of Fig. 2(c)] into the separate spectra of the two compounds, and the diffusion coefficients tell us their sizes.^63,64

Thus, IR-DOSY can be performed both in arrival and in depletion mode, and with similar accuracy. The advantage of the depletion mode becomes clear when one of the compounds has a diffusion coefficient that is so low that it becomes experimentally unfeasible to observe its arrival, as is the case with amyloids. In that case, IR-DOSY in depletion mode still makes it possible to obtain size-dependent species spectra.

We now demonstrate depletion-mode IR-DOSY on a mixture of BSA amyloids and monomers. The size of the BSA aggregates in our experiments ranges from nanometers to micrometers, leading to an extremely long time scale (≫1 day) for diffusion over the width of the sample-cell (4 mm). However, for the specific purpose of disentangling the spectral signatures of monomers and aggregates, it is sufficient to observe the time-dependent decrease of the monomer signal in the solution-filled half of the sample volume. The total time-dependent absorption in this part of the cell is the sum of that of the monomers and that of all the larger species (viz. the aggregates) in the sample that have diffusion constants that are so low that they are effectively immobilized. Thus, from a global analysis of the time- and frequency-dependent data, we can obtain the diffusion coefficient (and, hence, the size) of the monomers and the spectral signatures of both the monomers and the aggregates.

Figure 3 shows the FTIR spectra of BSA monomer (5–10 nm) and of a solution containing aggregates (∼200 nm), which was prepared by following the procedures described in the Methods section. The monomer spectrum consists of a broad band at ∼1655 cm⁻¹, with a shoulder at 1630 cm⁻¹, which are assigned in the literature to α-helices and unordered short segments between α-helixes.⁶⁵ The aggregate spectrum consists of a broad band at ∼1650 cm⁻¹, with a shoulder at ∼1615 cm⁻¹, which is due to the formation of inter-molecular β-sheets.⁶⁵ As was reported previously,⁶⁵ BSA aggregates show the typical band assigned to the beta-sheet structure, but the main spectral signature is a broad band around 1650 cm⁻¹, rather similar to the monomer spectrum, but slightly red-shifted and broader. This band is believed to be due to α-helical and unordered structures arising from partial denaturation.

The linear IR spectra of BSA amyloids and monomers are rather similar, and to enhance the contrast between their spectral signatures, we use second-derivative infrared spectra, in which the visibility of shoulders and peaks is enhanced compared to those of the conventional spectrum.^35,66 Figure 3(b) shows the second derivative (SD) of the infrared spectra shown in Fig. 3(a). In the monomer spectrum, we observe two well-separated minima at 1655 cm⁻¹ and at 1630 cm⁻¹, and a noisier minimum at 1678 cm⁻¹. In the aggregate spectrum, well-separated minima are present at 1615, 1640, 1650, and 1683 cm⁻¹. Although the second-derivative spectra have a lower signal-to-noise ratio, the weak signature assigned to the β-sheet in the normal IR spectrum of the aggregates shows up as a well-resolved minimum in the second-derivative spectrum.

We now study an ∼1:1 (w/w) mixed solution of BSA monomers and aggregates, using IR-DOSY, to disentangle the spectral signatures of the monomers and aggregates in the second-derivative spectrum. Figure 4 shows the time-dependent conventional IR spectrum recorded in depletion mode. Based on separate IR-DOSY experiments on the BSA monomer (see the supplementary material), this time range is sufficient to observe the diffusion of the monomeric form of the protein. In the conventional IR spectrum, we observe a decrease of the infrared absorption bands in the amide region, with the band at 1655 cm⁻¹ decreasing in intensity by ∼25%, as expected for complete monomer diffusion (and given the similar absorption intensities of the monomers and aggregates at this frequency). The time-dependent second derivative spectrum of the same experiment is shown in Fig. 4(b). Here, we observe a more significant change in intensity and a shift in the minimum at 1655 cm⁻¹, while the band at 1615 cm⁻¹, which is due to the beta-sheet structure of the aggregates, does not vary in intensity on the experimentally accessible time scale. In Fig. 4(c), we plot the second-derivative intensity vs time at the frequency of the main peak at 1655 cm⁻¹. By least-squares fitting the time- and frequency-dependent dataset with the solution of the diffusion equation (see the Methods section above), we obtain the diffusion coefficient of the monomer and the spectra associated with the monomer and the aggregates. In this fit, the diffusion constant of the aggregates was assumed to be negligible. The diffusion coefficient of the monomer obtained from the IR-DOSY experiment is similar to the one found using the IR-DOSY in arrival mode (see the supplementary material), as already mentioned, from the IR-DOSY experiment, but we can still obtain their spectral signature, which is shown in Fig. 4(d). Thus, we obtain a partial IR-DOSY spectrum, as shown in Fig. 5. The aggregate and monomer spectra obtained from the IR-DOSY experiment on the mixture agree well with those of the isolated species [Fig. 3(b)]: the monomer spectrum shows the same minima around 1655 and 1630 cm⁻¹, and weaker signature around 1610 and 1680 cm⁻¹, which have been assigned previously to aromatic side chains and turns in the monomer.⁶⁵ In the spectrum of the aggregates, we observe three minima at 1615, 1635, and 1650 cm⁻¹, in agreement with Fig. 3(b). Thus we see that IR-DOSY in depletion mode can separate the IR spectra of protein monomers and aggregates in a mixed solution in an unambiguous manner.

As was already mentioned in the Introduction, two-dimensional infrared spectroscopy has become a valuable lab tool in the study of amyloid formation and structure.^{11–13,15–19,36} The β-sheet diagonal and cross-peak features in the 2DIR spectrum of amyloids provide insight into their structure,^{11–13,15–19} and, in addition, excitonic effects in the 2DIR spectrum can give insight into the size of the amyloids.¹⁴ However, for mixed amyloid/monomer samples, 2DIR spectroscopy can become difficult due to the overlap of the amyloid and the monomer 2DIR spectra. This overlap problem is demonstrated for the case of BSA in Fig. 6, which shows the 2DIR spectrum (recorded in pump–probe mode) of a mixture of BSA amyloids and monomers, together with the 2DIR spectra of the pure monomer and amyloid. In the 2DIR spectrum of the mixture, the characteristic β-sheet cross-peak pattern³⁰ is largely hidden beneath the broad response of the monomers, in the center of the 2DIR spectrum. Just as in the case of the linear IR spectrum above, we can use the difference in size (and, hence, in the diffusion coefficient) of the amyloids and monomers to separate their contributions to the 2DIR spectrum.

We previously demonstrated such a 2DIR-DOSY spectroscopy in “arrival” mode for a mixture of small molecules.³⁸ Adding a third diffusion-coefficient axis to the 2DIR spectrum gives rise to three-dimensional spectra, in which 2DIR peaks of molecular species with different sizes show up as three-dimensional blobs (Fig. 5 of Ref. 38). In the present case, the diffusion constant of the monomer has already been obtained in the linear IR-DOSY experiments (see Sec. III B), so, we adopt a simpler strategy, shown schematically in the top part of Fig. 7. After injection of the sample solution and solvent, we wait for 24 h–a time that is sufficient for complete diffusion of the monomer over the sample volume to occur, but short enough that aggregate diffusion is still negligible. On the time scale of the experiment, the shift in the monomer/aggregate equilibrium due to the local changes in monomer concentration is negligible, because of the high stability of the aggregates and the slow formation rate of the aggregates (see supplementary material Fig. S4). The resulting distribution of monomers and amyloids at t = 24 h is shown schematically above Fig. 7(a): the monomer concentration is the same everywhere in the sample (viz., 50% of the initial concentration in the sample solution), whereas all the amyloids are still located in the bottom part of the sample volume (on the time scale of the experiment, the shift in the monomer/aggregate equilibrium due to the local changes in the monomer concentration is negligible, because of the high stability of the aggregates and the slow formation rate of the aggregates; see supplementary material Fig. S4). Hence, the 2DIR spectrum in the top part of the sample cell is the pure monomer 2DIR spectrum [Fig. 7(b)], and the 2DIR spectrum in the bottom part of the sample cell [Fig. 7(b)] is that of a mixture of amyloids and monomers, with exactly the same monomer concentration as in the top part [this monomer concentration is 50% reduced compared to the original sample solution, so, the spectra in Figs. 6(a) and 7(b) are not the same]. Hence, by subtracting the 2DIR spectra of Figs. 7(a) and 7(b), we obtain a different spectrum that is purely due to the amyloids, Fig. 7(c). This pure 2DIR spectrum of the amyloids obtained from the 2DIR-DOSY spectrum of the mixture is very similar to that of a solution of pure amyloid [Fig. 6(c)], and clearly shows the characteristic A_‖/A_⊥ beta-sheet cross-peak pattern,³⁰ which is difficult to discern in the raw 2DIR spectrum of the mixed sample solution [Fig. 6(a)].

An important advantage of the above method to obtain the pure amyloid 2DIR spectrum (compared to, e.g., measuring the monomer 2DIR spectrum in a separate sample and subtracting that from the mixture’s spectrum) is that it is not necessary to know the monomer or amyloid concentrations: after complete diffusion of the monomers, their concentration in the top and bottom parts of the sample cell is exactly the same, so in the subtraction of the 2DIR spectra in the top and bottom parts, the monomer contributions cancel out exactly. Furthermore, depending on the waiting time used, one can cancel out not only the contribution of the monomers to the 2DIR spectrum, but also that of oligomers: one only has to choose the waiting time sufficiently long that species to be canceled out have completely diffused over the sample volume, and the characteristic time for this process is L²/D, where L is the sample-channel width, and D = 6πηR is the diffusion constant, which is completely determined by the viscosity η of water and the size R of the molecules or oligomers. We also note that the waiting time can be reduced significantly by reducing the channel width: reducing the channel width from 4 mm to 400 μ will reduce the waiting time required by a factor of 100. Finally, combining IR-DOSY in depletion mode with the recently developed 2DIR method^67,68 to study proteins in water instead of D₂O should make it possible to also investigate amyloids in H₂O in the manner reported here.

We have shown that infrared diffusion-ordered spectroscopy in depletion mode provides a simple method to cleanly separate the monomer and amyloid contributions to the 1D-IR and 2D-IR spectra of mixed monomer/amyloid samples. The depletion mode is necessary because the diffusion of the large aggregates takes place on a time scale that is experimentally not accessible. The only drawback is that only the size of monomers, and not that of the aggregates, is obtained in these IR-DOSY experiments. However, if the particles under study are so large that this becomes a problem, their size can generally be determined relatively easily with other methods such as dynamic light scattering or microscopic imaging.

IR-DOSY should also be applicable if, in addition to amyloids and monomers, there are oligomers present in the sample solution. In this case, the temporal evolution of the (2D)IR spectrum cannot be described by the diffusion of a single species (as was the case in the BSA monomer/amyloid mixture studied here), but requires global fitting with two or more diffusion coefficients and their associated (2D)IR spectra, with the diffusion coefficients providing information on the size of the oligomers and the associated (2D)IR spectrum on their structure. We have demonstrated previously that such multi-component samples can be analyzed with IR-DOSY,³⁸ and we are working on more efficient algorithms to determine the number of species directly from the multidimensional dataset and to perform the least-squares fitting. Such an analysis, in terms of a discrete set of sizes and spectra, might not work in the case of more complex mixtures containing monomers, oligomers with varying sizes, and medium to large amyloids, which requires an analysis in terms of continuous size (and hence diffusion-coefficient) distributions. We are currently developing the algorithms for analyzing the (2D)IR-DOSY data of such samples using the sophisticated analytical methods developed for NMR-DOSY (such as the constrained-regularization^40,69,70 and the maximum-entropy⁷¹ methods), which have proven successful in characterizing samples with continuous size distributions. We believe that, after such further development, 2DIR-DOSY can become a valuable technique to reveal the presence and structure of amyloidic oligomers, and, so help improve our understanding of the formation and structure of amyloids.

The supplementary material consists of: conventional IR spectra of acetone and dialanine; IR-DOSY of BSA monomer in heavy water; details of the fitting procedures.

We thank Michiel Hilbers, Wim Roeterdink, and Hans Sanders for technical support. F.C. acknowledges financial support from The Netherlands Organization for Scientific Research (NWO) (Grant No. 680-91-13). This publication is part of the project (Project No. VI.Veni.212.240) of the research program NWO Talent Program Veni 2021, which is financed by the Dutch Research Council (NWO).

The authors have no conflicts to disclose.

Giulia Giubertoni: Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Federico Caporaletti: Data curation (supporting); Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting). Rianne Van Diest: Formal analysis (supporting); Investigation (equal); Visualization (equal); Writing – original draft (equal). Sander Woutersen: Conceptualization (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.