Heat flows enrich prebiotic building blocks and enhance their reactivity

List of substances used

For a list of abbreviations, see Supplementary Table 1. Nucleobases A and T, nucleosides Ado, Guo and Urd, nucleotides (2′-AMP, 3′-AMP, 5′-AMP, 2′,3′-AMP, 3′,5′-AMP, 5′-CMP, 2′,3′-CMP, 5′-GMP, 3′,5′-GMP and 5′-UMP), deoxyribonucleotides (dAMP, dCMP, dGMP and dTMP), 2-aminoazoles (2AO, 2AI and 2AT), oligomers of Gly (Gly, GlyGly and GlyGlyGly) and TMP were purchased from Sigma-Aldrich (USA). Nucleobase C was purchased from Carl Roth GmbH (Germany) and nucleobases U and G and nucleosides L-Cyd and D-Cyd were purchased from Biosynth Carbosynth (UK). Nucleotides 2′-CMP, 3′-CMP, 3′,5′-CMP, 2′,3′-GMP, 2′,3′-UMP and 3′,5′-UMP were purchased from Biolog (USA). For the proteogenic amino acids, an Amino Acid Mixture was used (L4461, Promega, USA) and compared with individual amino acids (Sigma-Aldrich, USA). For experiments on non-proteogenic amino acids, a mixture of amino acids and small molecules was used (A9906, Sigma-Aldrich, USA).

Experimental procedure: preparation of heat flow cells

Bold numbers refer to the corresponding encircled numbers in Extended Data Fig. 1. The FEP foils (5; Holscot, the Netherlands) defining the microfluidic structure were cut by an industrial plotter device (CE6000-40 Plus, Graphtec, Germany). The cutout was then sandwiched between two sapphires (4 and 6; KYBURZ Switzerland, Switzerland) with thicknesses of 500 µm (cooled sapphire) and 2,000 µm (heated sapphire), respectively. The sapphires were previously coated with a hydrophobic coating (ProSurf MT-5, Surfactis, France) to avoid interaction with the sample and facilitate sample extraction. The cooled sapphire (4) has four laser-cut holes with a diameter of 1 mm. The sapphire–FEP–sapphire block was then placed on an aluminium base (2), covered by a heat-conducting foil (3; graphite, 25 µm, 1,600 W mK⁻¹, Panasonic, Japan) and held in place by a steel frame (7), which is connected to the aluminium base by six torque-controlled steel screws for homogeneous force distribution.

The height of the chamber was measured with a confocal micrometer (CL-3000 series with CL-P015, Keyence) at three positions (bottom, middle and top) to ensure a homogeneous thickness of the chamber. The average of these three measurements is later used in the numerical model to determine the Soret coefficient for each experiment. An Ohmic heating element (9) was mounted on the steel frame with torque-controlled screws and connected to the heated sapphire with another heat-conducting foil (8; EYGS0811ZLGH, graphite, 200 µm, 400 W mK⁻¹, Panasonic).

For the microfluidic connections, we used (Techlab, Germany): connectors (UP P-702-01), end caps (UP P-755), screws (VBM 100.823-100.828), ferrules (VBM 100.632) and tubings (Teflon (FEP), KAP 100.969). As syringes, we used (ILS, Germany, bought from Göhler-HPLC Syringes, Germany) 2606814 and 2606035. Chambers were pre-flushed using low-viscosity, fluorinated oil (3M Novec 7500 Engineered Fluid, USA) to check for tightness and push out residual gas inclusions.

For the 100 mM phosphate buffer used in Extended Data Fig. 5c, we dissolved 584 mg NaH₂PO₄·H₂O and 819 mg Na₂HPO₄ in 100 ml water, resulting in a solution of pH 7. Each mixture was chosen to cover a broad range of molecules while still producing separable peaks in the HPLC measurements (see Extended Data Fig. 1 and Supplementary Fig. 14).

The sample (42.5 µl) was loaded into the chamber with the help of a syringe filled with fluorinated oil to avoid the inclusion of air bubbles. After loading the sample, the tubings were closed with end caps.

The chamber was then mounted to a cooled aluminium block connected to a cryostat (TXF-200 R5, Grant, UK). The heaters were connected to a 400-W, 24-V power supply controlled by means of Arduino boards with a customized version of the open-source firmware Repetier, initially designed for 3D printing. The cryostat and the heaters were set to the desired temperatures. Temperatures were measured on the sapphires with a heat imaging camera (ShotPRO, Seek Thermal, USA), giving cold and hot temperatures per experiment.

All elements were optimized using finite element simulation to provide a homogeneous temperature profile.

To stop the experiment, heaters and the cryostat were turned off and the chamber was frozen at −80 °C for at least 15 min. This allowed us to open the chamber and cut the frozen interior into four fractions of equal volume. Fewer fractions would lower resolution by averaging over a larger fraction of the chamber (Extended Data Figs. 1 and 2), whereas more fractions would suffer a larger positional error from manual cutting (for a comparison of 12 versus four fractions, see Extended Data Fig. 2).

For experimental network experiments, a modified cutout was used (shown in Fig. 4). Before mixing with concentrated amino acid standard, the water used for dilution was degassed by heated stirring under vacuum. The sandwich was filled with sample solution before assembly to minimize air inclusions and dilution effects from remaining water. To start the experiment, inlets and outlets were connected to tubings (250 µm inner diameter to reduce air inclusions, KAP 100.966, Techlab, Germany) filled with sample solution and connected to syringes driven by syringe pumps (Nemesys, CETONI, Germany). As soon as the temperature gradient (16 K between 30 °C and 46 °C and 10 K between 32 °C and 42 °C) was established, the inlet was supplied with a continuous flow of 1 nl s⁻¹ and the outlets with 0.5 nl s⁻¹ each. For the recovery of fractions of the three individual chambers, a slightly modified extraction procedure was used. Because the time required to extract 12 fractions (four fractions per chamber) is increased compared with the four fractions used above, we added dry ice around the frozen block to reduce dilution by condensation. Also, we directly pipetted the recovered fraction into 7 eq. of borate buffer for pre-column derivatization.

pH measurements

For pH measurement and adjustment, we used a Versa Star Pro device equipped with an 8220BNWP micro pH electrode (Thermo Fisher Scientific, USA). The fundamental effect of heat-flow-driven pH gradients was shown previously⁴¹. However, given the low concentrations of solutes used in this work, the obtained pH differences remained weak (delta pH between 0 and 1 units).

LC + MS measurements

Two different LC systems were used for the various analysis methods for all molecules of interest. Also, for glycine dimerization and detection of non-proteogenic amino acids, an orbitrap MS was used. A method was created and optimized for each set of species to ensure the best separation. In the following, these will be listed for the molecules involved.

System 1: the system consists of a Vanquish Flex (VF-S01-A), a binary pump (VF-P10-A-01), a heated column compartment (VH-C10-A) and a variable-wavelength detector (VF-D40-A; all Thermo Fisher Scientific, USA). As column, we used Symmetry C18 (3.5 µm pore size, 2.1 mm diameter, 150 mm length, 100 Å particle size, WAT106005) (Waters, USA). As eluents: eluent A: LC-H₂O (0.1% v/v formic acid); eluent B: LC-acetonitrile (0.1% v/v formic acid). For all methods, we used a flow of 0.3 ml min⁻¹, a column temperature of 30 °C (still air) and detection of the UV absorption at 260 nm with 50 Hz. All methods are followed by a washing step at 40% B for 1 min and equilibration at starting concentration for 6 min.

Detailed protocols for the respective mixtures of compounds:

• Nucleobases: isocratic elution for 3 min at 0% B, then increase to 10% B over 2 min.

• Nucleosides: isocratic elution for 3 min at 0% B, then increase to 5% B over 4 min.

• Cytidine monophosphates: isocratic elution for 5 min at 0% B, then increase to 5% B over 2 min.

• Adenosine monophosphates: start with 0% B, then increase to 5% B over 12 min.

• 5′-ribonucleotides: isocratic elution for 10 min at 0% B.

• 2′,3′-cyclic ribonucleotides: start with 0% B, then increase to 4.5% B over 15 min.

• 3′,5′-cyclic ribonucleotides: start with 0% B, then increase to 6.8% B over 9.5 min.

• Deoxyribonucleotides: start with 0% B, then increase to 3% B over 8 min, then increase to 5% B over 4 min and 7.5% B over 1 min.

• Dimerization of glycine (variation of TMP): pre-column derivatization as described below. Isocratic elution for 10 min at 4% B.

MS-based methods: we carried out mass analysis using a Q Exactive Plus Orbitrap HR/AM (Thermo Fisher Scientific, USA), using positive ionization with a resolution of 70k, an AQC target of 3 × 10⁶ and a maximum IT of 200 ms. On the HESI source, a sheath gas flow rate of 2 was set, a spray voltage of 2.9 kV, 320 °C capillary temperature, 50 °C auxiliary gas heater temperature and a S-lens RF level of 50. For analysis, the main isotope mass ±0.075 m/z was extracted. For LC-MS methods, we used the same settings for the LC as above.

Detailed protocols for the respective mixtures of compounds:

• Dimerization of glycine (time variation and heat flow chambers): isocratic elution for 5 min at 4% B.

• Mixture of non-proteogenic amino acids, proteogenic amino acids, dipeptides and other molecules: isocratic elution for 5 min at 4% B.

System 2: the system consists of a Vanquish Core (VC-S01-A-02), a quaternary pump (VC-P20-A-01), a heated column compartment (VC-C10-A-03), a diode array detector (VC-D11-A-01) and a fluorescence detector (VC-D50-A-01; all Thermo Fisher Scientific, USA).

Separation of 2AO, 2AI and 2-aminothiazole (2AT): following a method from the literature⁵⁰, we applied isocratic elution for 5 min with 90% LC-H₂O, 10 mM ammonium formate and 10% LC-acetonitrile. The flow rate of 1.5 ml min⁻¹ was applied on the column InertSustain Amide (100 Å, 5 µm, 4.6 × 150 mm, GL5020-88631) (GL Sciences, Japan), with a column temperature of 40 °C and UV detection at 225 nm.

Separation of all 20 proteogenic amino acids: for the separation of amino acids, we adapted a method from the literature⁵¹. Following that protocol, we first pre-column derivatized our sample. This was done by mixing 14 µl of 50 mM borate buffer (28341, Thermo Fisher Scientific, USA) with pH adjusted to 8.8 and 2 µl of our sample. Then, 4 µl of freshly prepared 6-aminoquinolyl-N-hydroxysuccinimidyl carbamate (AQC, S041, Synchem, USA) and 4 mg ml⁻¹ in anhydrous acetonitrile (43166, Alfa Aesar, USA) were added and mixed with a pipette. The samples were then incubated for 10 min at 55 °C and stored in the autosampler at 5 °C before injection. We used the column ACCLAIM Vanquish C18 (2.2 µm, 2.1 mm × 150 mm) (Thermo Fisher Scientific, USA). Eluent A was LC-H₂O (W6-212, Fisher Scientific, USA) + 50 mM ammonium formate (17843, Honeywell, USA) + 0.8% v/v formic acid (A117-50, Fisher Scientific, USA) and eluent B was LC-acetonitrile (A955-212, Fisher Scientific, USA). We set the column temperature at 45 °C (still air) and simultaneously acquired the UV absorption signal at 260 nm and the fluorescence signal (FLD) under excitation at 266 nm and emission at 473 nm (both with 50 Hz). Elution was done with a flow of 0.65 ml min⁻¹. Best separation was achieved using the following elution protocol: isocratic elution with 0.5% B between 0 and 0.548 min, increase to 5.2% B until 3 min, to 9.2% B until 8.077 min, to 14% B until 8.626 min, maintain at 14% B until 9.5 min, increase to 19.2% B until 11.227 min, to 19.5% B until 13.696 min, to 90% B until 14.4 min and finally lower to 0.5% B and equilibrate for 6 min. For peak identification, samples of the individual amino acids were prepared according to the same protocol. Calibrations with different concentrations of amino acid standard solutions were done using the same volume as for experimental samples and for each set of measurements. After a comparison of the results of these calibrations in UV and FLD measurements, we found that the two detection methods worked equally. Because background levels in FLD are very low, we chose to proceed with the FLD channel except for tryptophan, which does not fluoresce, and tyrosine, which is known to produce intermediary products⁵¹. For Fig. 2 and Extended Data Fig. 3, we injected each sample twice and used the average of these two injections.

Treatment of data

Integral values of all experiments were analysed using Python. The accumulation profiles (Fig. 2a and Extended Data Figs. 2a and 8) over the chamber height were internally normalized for each species (see equation (1)) to compare enrichments despite the different initial concentrations owing to the use of different mixtures. Inherently, this makes the system compatible with mixtures with different concentrations of species. The normalized concentration [A]_j,k for each individual measurement was calculated for species A for fraction j ∈ {1, 2, 3, 4} and replicate number k ∈ {1, 2, 3}.

with the HPLC-measured concentration [A]_j,k,HPLC obtained from the integral values presented in Supplementary Tables 5–65. For the detection of the dimerization of glycine, we used a quadratic function of \((\sqrt{{\rm{response}}})\) to account for the nonlinear response (Supplementary Information pages 166–167). For all other species, we used an external linear calibration (Supplementary Tables 3 and 4 and examples on pages 135–165 of the Supplementary Information). The linear calibration intrinsically emphasizes that differences in calibration or absolute concentration do not change the resulting enrichments.

To compare two species inside a pool of molecules (as shown in the heat maps in Figs. 2 and 3 and Extended Data Figs. 2–5), we calculated the concentration ratio of species A compared with species B in the top (j = 1) and respective bottom (j = 4) fraction for each combination of species and per individual experiment. We then calculated the average of the triplicate experiments.

The averaging over the triplicates is done only after the calculation of the concentration ratios. This is necessary because the temperature gradients between the replicate experiments differ slightly (1–2 K). This affects the concentrations of all species present in the respective mixture of a replicate equally (Supplementary Fig. 9), so that a calculation of the concentration ratio only after averaging the species concentrations would lead to a distortion of the enrichment value actually present in the heat flow chamber. The error maps shown in the Extended Data Figures and Supplementary Figures were determined by calculating the s.d. of the enrichments [A]_j,k/[B]_j,k of all repeats from the average value described in equation (2).

Control experiments were done in parallel to check for possible degradation of individual compounds. For this purpose, one solution each was incubated at the lower or higher temperature occurring in the heat flow chamber in bulk for the duration of the accumulation experiments, but no substantial selective degradation was observed.

For the enrichment of an individual species compared with a pool of molecules (Supplementary Table 2), we normalized the individual experiments as explained previously. Then, we calculated the average concentration \({\bar{c}}_{j,k}\) of the pool of all species k (per experiment) in the top (j = top) and respective bottom (j = bot) fraction for species i ∈ {1,…, S}.

By comparing the concentration of the species with the mean concentration of the pool \({\bar{c}}_{j}\), we were able to determine the enrichment of species A in fraction j (per experiment), which we averaged over the triplicate measurements (Supplementary Table 2).

Determination of the thermophoretic strength (Soret coefficient) of prebiotic substances in complex mixtures

A new measurement approach was necessary to determine the thermophoretic properties of the very small and highly diluted prebiotically relevant chemicals mixed together in networks of interconnected rock fractures. Previous methods require either fluorescent labelling of the measured substance, which distorts the thermophoretic properties of small molecules, or high substance concentrations, which—in our example—would strongly change the pH value and, thus, would not be representative for the prebiotic context of diluted solutions. At the same time, only one substance can be measured at a time. To overcome these difficulties, we use the same setup that mimics the geological scenario of heat flows through thin rock fractures that has also been used in previous studies^40,42 (Extended Data Fig. 1). The microfluidic structure is prepared as explained earlier. The combination of thermal convection and thermophoresis allows the solutes to be separated along its height of 50 mm much more effectively than would be possible by thermophoresis alone⁵². This spatial separation allows the entire chamber to be frozen. The content of the chamber is then cut into four equal-sized pieces (12.5 mm height each) using a scalpel. All pieces are analysed separately by HPLC (Extended Data Fig. 1), resulting in an average concentration \({\overline{[{\rm{A}}]}}_{j}=1/3\sum _{k}{[{\rm{A}}]}_{j,k}\) of species A at position j, with [A]_j,k as introduced in equation (1).

We assumed diffusive mobilities of D_i,AA ≈ 800 µm² s⁻¹ for amino acids⁵³ and D_i,AA ≈ 1,400 µm² s⁻¹ for the 2-aminoazoles and nucleotide components⁵⁴. This approximation is reasonable because the experiment is close enough to steady state after 18 h and, hence, the fitted thermophoretic strength S_T,i depends only marginally on D_i (Supplementary Fig. 8).

To determine the thermophoretic strength from these datasets, we first create a 2D model of the same height (50 mm) and thickness (0.17 mm) as in the experiment using finite element methods (COMSOL 5.4). In this model, the temperatures of the sidewalls are set according to the experiment. To determine the thermal convection of the solvent, we solve the Navier–Stokes equation

and the continuity equation in the steady-state case

in which ρ denotes solvent density, u the solvent velocity vector, I the unit vector, p the pressure, µ the dynamic viscosity and g the gravitational acceleration. The results show a laminar convection flow u inside the cell, which is coupled to the solute movement as it drags it with it.

This is achieved by solving the time-dependent drift–diffusion equation

with the local solute concentration c_i,num of species i (initial concentration c_i,num,0 = 1), its diffusive mobility D_i and Soret coefficient \({S}_{T,i}\equiv \frac{{D}_{T,i}}{{D}_{i}}\), which is defined using its thermophoretic mobility \({D}_{T,i}=-\frac{{v}_{T,i}}{\nabla T}\), including the thermophoretic drift velocity v_T,i. D_i is determined approximately by literature values after making sure that it does not change the determination of S_T,i (Supplementary Fig. 8). To obtain S_T,i for each species of the mixture, we fitted the numerical results of the averaged concentrations of all four volume fractions V_j, \({\bar{c}}_{i,j,{\rm{num}}}={\int }_{{V}_{j}}{c}_{i,{\rm{num}}}{\rm{d}}V\) to the average concentrations obtained by HPLC \({\overline{[{\rm{A}}]}}_{j}\), i = A. Although external control of COMSOL is possible to include it directly in the fitting algorithm, we chose to first solve equations (5)–(7) for a wide range of parameters and then use this dataset to linearly interpolate the experimental results, yielding the respective Soret coefficient of the solute much faster and with good precision. The parameter range covered different temperature gradients (ΔT [K]: 0, 5, 10, 20, 25), thermophoretic mobilities \(({D}_{T}\left[{\times 10}^{-12}\,{{\rm{m}}}^{2}\,{{\rm{sK}}}^{-1}\right]:1,\,2.5,\,5,\,10,\,20,\,40,\,60,\,80)\) and diffusive mobilities \((D\left[\times {10}^{-12}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}\right]:1,700,\,1,400,\,1,100,\,700)\) for 100 equidistant time points between 0 h and 24 h, which results in 192 different concentration profiles at each time point (see model file SimpleSim_2022_03_08_new.mph, resulting in the data file 2022_03_08_2DSim.dat). Using a custom-made LabVIEW program incorporating Levenberg–Marquardt algorithms (see SingleNTD_TrapFitter_V1-6.llb), we varied the D_T,i with fixed, experimentally obtained values for ΔT, D_i until an optimal fit between the numerical \(({\bar{c}}_{i,j,{\rm{num}}})\) and experimental \(({\overline{[{\rm{A}}]}}_{i})\) concentration profiles was found for each species. This procedure was repeated for all experimental repeats (triplets), after which the obtained values for S_T,i were averaged and a s.d. was obtained (Extended Data Table 1). Steady-state concentration profiles shown in Extended Data Fig. 2c,d and Supplementary Fig. 6 were obtained by calculating

in which β denotes the volume expansion coefficient of water, α the distance between the hot and the cold sides of the heat flow chamber, η the dynamic viscosity of water and y the space coordinate along the height of the chamber.

Modelling of a system of connected cracks/heat flow chambers

To determine the concentration profiles in a system of interconnected heat flow chambers from the previously determined Soret coefficients, we first had to calculate the behaviour of all species in a single heat flow chamber under various boundary conditions, such as temperature difference and flow rates.

For this, we created a 3D chamber in COMSOL with a height of 200 mm, a width of 60 mm and a thickness of 0.17 mm. The chamber has an inflow channel at its top end and an outflow channel at its top and bottom ends (see supplied COMSOL file SingleNTD_v7_simpleGeo_Large_allSTs_60mmWide200mmHighTrap.mph, yielding data file SingleNTD_v7_simpleGeo_Large_60x200mmWide.dat). We then solve equations (5)–(7) in the steady-state case, calculating first the laminar flow and then the concentration distribution of the solute. Solutions are generated for all combinations of boundary conditions, that is, for different temperature differences (ΔT [K]: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), diffusive mobilities \(({D}_{i}\left[\times {10}^{-12}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}\right]:\,800,\,1,400)\), Soret coefficients \(({S}_{T,i}\left[\times {10}^{-3}\,{{\rm{K}}}^{-1}\right]:\,1,\,4,\,8,\,12)\) covering the range given by Extended Data Table 1, inflow volume rates \(({Q}_{{\rm{in}}}[\times {10}^{-2}\,{\rm{nl}}\,{{\rm{s}}}^{-1}]:\,1,\,5,\) \(10,\,50,\,100,\,500,\,1,000)\) and outflow rates through the bottom outflow channels as ratios of the inflow volume rate (Q_out,bot [%]: 1, 5, 10, 20, 50, 100). Because the random distribution of channels in the modelled network of heat flow chambers results in channels with flow rates smaller and larger than the assumed 1 nl s⁻¹ (Supplementary Fig. 10), we appropriately set the range of simulated flow rates Q_in. The outflow rate of the upper channel is determined from the mass conservation of the incompressible aqueous solution. This parametric sweep yields a dataset of 7,332 concentration profiles.

In the second step, we use this extensive dataset to extrapolate the behaviour of interconnected heat flow chambers. For this, we set up the structure of the chamber system using a self-made LabVIEW program that fills a 3D array (see NetworkSimulation_v3-1.llb with the settings shown in GridSimSettings.png and the program state saved in MCA31_2023_07_24.dat), representing Z systems of a 2D matrix of chambers. The program assigns one inflow and a maximum of two outflow channels to each of the chambers in a random fashion. Each chamber can only connect to chambers inside the same row or one row above or below, simulating a realistic system. The inflow and outflow rates are set up according to the mass conservation of the solvent. After setting the temperature difference ΔT and defining all species i of a mixture by setting D_i and S_T,i, every chamber of the system has a fully defined set of boundary conditions that can be retrieved from the previously calculated dataset. Solutions of parameter values (ΔT, D_i, S_T,i, Q_in, Q_out,bot) that are set between calculated parameters are linearly interpolated from the above dataset. The concentration profiles of each chamber c_m,n,i at matrix position m and n are then renormalized column by column with the total inlet concentration \({c}_{{\rm{in}},m,n,i}=\frac{{\sum }_{l}{c}_{{\rm{in}},m,n,i,l}\times {Q}_{{\rm{in}},l}}{{Q}_{{\rm{in}}}}\), with the concentrations c_in,m,n,i,l and volume rates Q_in,l (Q_in = Σ_l Q_in,l) of each of the l_max inflow channels l. The outflow concentrations c_out,m,n,i of each species are calculated by the finite element simulations and assigned to the inlet concentrations c_in,m,n+1 of the next column. The renormalized averaged concentrations c_bot,m,n,i of the bottommost 1.5 mm volume fractions are then plotted in Fig. 4e–i. The maximum pairwise enrichments shown in Fig. 4g–i are calculated by taking the median of the ten maximum concentration ratios of two species A and B: \({\rm{median}}\,\left(\mathop{\max }\limits_{\forall m,n}\left(\frac{{c}_{{\rm{bot}},m,n,i={\rm{A}}}}{{c}_{{\rm{bot}},n,m,i={\rm{B}}}}\right)\right).\) The median is used to avoid overestimating enrichment ratios by numerical noise and outliers.

Owing to the routing of the connecting channels according to the above rules as well as the random error of the measured Soret coefficients, the numerical modelling of the heat flow chamber network also has associated errors. To quantify these, we repeated the respective simulations (that is, N_x by N_y chambers for Z systems) at least three times and calculated the mean and corresponding s.d. from the resulting enrichment values as shown in Extended Data Fig. 6. To determine the random error for the Soret coefficients shown in Extended Data Table 1, we determined the correlation of the Soret coefficients for all species present in the respective mixtures between measurement replicates (Supplementary Fig. 9). The systematic error is mainly caused by small differences in the temperature gradient between the individual replicates and equally affects the Soret coefficients of all species contained in a mixture. By fitting this correlation linearly, this systematic error can thus be determined by the deviation of the slope of the linear fit to the slope 1. Accordingly, the s.d. of the fit yields the random error (for example, because of the integral determination of the HPLC peaks, fluctuations of the column temperature during the HPLC run etc.). In the repeat runs of the network modelling for the error determination, the Soret coefficients for the species involved were, therefore, randomly chosen according to a Gaussian distribution around the average value shown in Extended Data Table 1 with a s.d. corresponding to the random error determined above. Thus, the error from the Soret coefficients could be numerically propagated and shown in the error matrices in Extended Data Fig. 6.

Glycine-dimerization experiments

For the experiments shown in Fig. 5, we closely followed the protocol for the dimerization of glycine¹¹.

We prepared triplicates of 200 µl of a solution containing 10 mM and 100 mM glycine and various concentrations of TMP, ranging more than three orders of magnitude (0.10, 0.25, 0.63, 1.6, 4.0, 10, 25, 63, 158 mM). The mixture was then adjusted to pH 10.5 using NaOH and heated at 90 °C for 16 h. Using the LC protocol described above, we then analysed the yield of dimerization by comparison against standards.

The same was done for a mixture of 100 mM glycine and 100 mM TMP with measurements at different time points between 0 and 120 h.

Modelling of reactions for network extrapolation

To model the dimerization of glycine driven by TMP in large networks of connected heat flow chambers, we first determined all required rate constants k₁–k₅ of the simplified reaction scheme^11,12:

Here the single glycine in equation (9) is first phosphorylated (activated) by TMP. The activated glycine can then react with another glycine in equation (11) and form the GlyGly dimer. Both the activated glycine and the glycine dimer can degrade to glycine by hydrolysis (equations (10) and (13)). For TMP, a possible hydrolysis was also taken into account (equation (12)). For the determination of k₁–k₅, we measured the product concentrations for a range of initial concentrations ([Gly]_init = 10 mM, 100 mM at [TMP]_init = 100 mM, T = 90 °C, initial pH 10.5) at different time points (Fig. 5c and Supplementary Table 17). We determined the rate k₄ separately with a solution of only [TMP]_init = 0.2 mM in water under the same temperature and pH conditions to allow a better stability of the numerical fit described below for the remaining reaction constants k₁–k₃, k₅ and measured the decrease in TMP concentration over time (Supplementary Fig. 13 and Supplementary Table 19).

For the determination of reaction rates, either the full reaction system (equations (9)–(13)) or only the TMP hydrolysis (equation (12)) was implemented as a differential equation system in a 0D model (COMSOL 5.4, filename: MultiDFitGlyRctnModel_20230714_Modellv1.mph, yielding the data file 2023_07_14_MultiDFit_modelV1.dat, analysed with FreeFitter_V4-6.llb):

The solution of this reaction system was then solved for the initial concentrations [Gly]_init = 100 mM and 10 mM and [TMP]_init = 0.0001 M, 0.000251189 M, 0.000630957 M, 0.00158489 M, 0.00398107 M, 0.01 M, 0.0251189 M, 0.0630957 M and 0.158489 M (corresponding to the experiment in Fig. 5c) under variation of the rate constants (Supplementary Table 18). The solution concentrations were determined for the reaction times 0–122 h. For the separate determination of k₄ mentioned above, the initial concentration [TMP]_init = 0.2 mM was chosen according to the experiment (Supplementary Fig. 13). With the extensive datasets obtained in this way, we were able to simultaneously fit all of the experimental data points to a common parameter set with a custom-made LabVIEW program using a Levenberg–Marquardt algorithm (see FreeFitter_V4-6.llb). The rates obtained in this way are: k₄ = 3.5 × 10⁻⁷ s⁻¹ (±19%) and k₁ = 1.0 × 10⁻³ Ms⁻¹ (±41%), k₂ = 1.1 × 10⁻⁴ s⁻¹ (±10%), k₃ = 9.7 × 10⁻⁵ Ms⁻¹ (±21%) and k₅ = 1.2 × 10⁻⁶ s⁻¹ (±20%). Although the error bars in Fig. 5c indicate the s.d. from triple replicate experiments, we have determined the error of the modelling described above by a stochastic approach. We calculated the resulting product concentrations 500 times for the fitted set of rate parameters, choosing the rates from the fit using a weighted random number generator with a Gaussian probability distribution with the s.d. according to the previously given rate error. The shaded area contains 68.27% of all these 500 runs, equivalent to one s.d. To determine the reaction yield in the network of interconnected heat flow chambers, we proceeded as described next.

Using the rate constants thus obtained, we calculated another dataset using COMSOL as the result of the equations (14)–(17), in which we varied the initial concentrations [Gly]_init and [TMP]_init over a range of 1 × 10⁻¹⁰ M to 1 M each, with five concentrations per decade (see GridReactionGlyRctnModel_20230717_Modellv1.mph, resulting in data file 2023_07_17_ModelV1ForGridSimFig5.dat for analysis with FreeFitter_V4-6.llb). The product concentrations were determined over a reaction time of 120 h (divided into 20 time points). We were thus able to map the TMP and Gly concentrations obtained from network modelling described above (using the Soret coefficients from Extended Data Table 1) in each individual heat flux chamber to this reaction dataset and thus determine the amount of product obtained and the reaction yield (Fig. 5d–f). The errors given in Fig. 5e,f were calculated stochastically as described above by recalculating the network model ten times, choosing the Soret coefficients with a weighted random generator with Gaussian probability distribution at one s.d. of the previously determined random error.

The network continuously provides the reactants according to the accumulation characteristics shown in Fig. 4, taking into account the high Soret coefficient of TMP S_T,TMP ≈ 7 × 10⁻³ K⁻¹. As shown in Fig. 5c, right, the maximum reaction yield is already reached after 16 h, but the relaxation time of a chamber with throughflow in the range Q_flow = 1–10 nl s⁻¹ is expected to be on the order of Volume_chamber/Q_flow ≈ 10¹–10² h. Therefore, to sufficiently account for the hydrolysis included in the model, we calculate the product yield after 120 h reaction time. The reactions are each assumed to occur in a reaction volume at the bottom of the respective chamber in which the reactants are most concentrated (Supplementary Fig. 11a).