[ad_1]
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), 2aminoazoles (2AO, 2AI and 2AT), oligomers of Gly (Gly, GlyGly and GlyGlyGly) and TMP were purchased from SigmaAldrich (USA). Nucleobase C was purchased from Carl Roth GmbH (Germany) and nucleobases U and G and nucleosides LCyd and DCyd 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 (SigmaAldrich, USA). For experiments on nonproteogenic amino acids, a mixture of amino acids and small molecules was used (A9906, SigmaAldrich, 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 (CE600040 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 MT5, Surfactis, France) to avoid interaction with the sample and facilitate sample extraction. The cooled sapphire (4) has four lasercut holes with a diameter of 1 mm. The sapphire–FEP–sapphire block was then placed on an aluminium base (2), covered by a heatconducting foil (3; graphite, 25 µm, 1,600 W mK^{−1}, Panasonic, Japan) and held in place by a steel frame (7), which is connected to the aluminium base by six torquecontrolled steel screws for homogeneous force distribution.
The height of the chamber was measured with a confocal micrometer (CL3000 series with CLP015, 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 torquecontrolled screws and connected to the heated sapphire with another heatconducting foil (8; EYGS0811ZLGH, graphite, 200 µm, 400 W mK^{−1}, Panasonic).
For the microfluidic connections, we used (Techlab, Germany): connectors (UP P70201), end caps (UP P755), screws (VBM 100.823100.828), ferrules (VBM 100.632) and tubings (Teflon (FEP), KAP 100.969). As syringes, we used (ILS, Germany, bought from GöhlerHPLC Syringes, Germany) 2606814 and 2606035. Chambers were preflushed using lowviscosity, 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_{2}PO_{4}·H_{2}O and 819 mg Na_{2}HPO_{4} 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 (TXF200 R5, Grant, UK). The heaters were connected to a 400W, 24V power supply controlled by means of Arduino boards with a customized version of the opensource 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^{−1} and the outlets with 0.5 nl s^{−1} 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 precolumn 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 heatflowdriven pH gradients was shown previously^{41}. 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 nonproteogenic 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 (VFS01A), a binary pump (VFP10A01), a heated column compartment (VHC10A) and a variablewavelength detector (VFD40A; 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: LCH_{2}O (0.1% v/v formic acid); eluent B: LCacetonitrile (0.1% v/v formic acid). For all methods, we used a flow of 0.3 ml min^{−1}, 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): precolumn derivatization as described below. Isocratic elution for 10 min at 4% B.
MSbased 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^{6} 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 Slens RF level of 50. For analysis, the main isotope mass ±0.075 m/z was extracted. For LCMS 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 nonproteogenic 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 (VCS01A02), a quaternary pump (VCP20A01), a heated column compartment (VCC10A03), a diode array detector (VCD11A01) and a fluorescence detector (VCD50A01; all Thermo Fisher Scientific, USA).
Separation of 2AO, 2AI and 2aminothiazole (2AT): following a method from the literature^{50}, we applied isocratic elution for 5 min with 90% LCH_{2}O, 10 mM ammonium formate and 10% LCacetonitrile. The flow rate of 1.5 ml min^{−1} was applied on the column InertSustain Amide (100 Å, 5 µm, 4.6 × 150 mm, GL502088631) (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^{51}. Following that protocol, we first precolumn 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 6aminoquinolylNhydroxysuccinimidyl carbamate (AQC, S041, Synchem, USA) and 4 mg ml^{−1} 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 LCH_{2}O (W6212, Fisher Scientific, USA) + 50 mM ammonium formate (17843, Honeywell, USA) + 0.8% v/v formic acid (A11750, Fisher Scientific, USA) and eluent B was LCacetonitrile (A955212, 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^{−1}. 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^{51}. 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}.
$${[{\rm{A}}]}_{j,k}=\frac{{[{\rm{A}}]}_{j,k,{\rm{HPLC}}}}{\frac{1}{4}{\sum }_{j}{[{\rm{A}}]}_{j,k,{\rm{HPLC}}}}$$
(1)
with the HPLCmeasured 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.
$$\overline{{[{\rm{A}}]}_{j}/{[{\rm{B}}]}_{j}}1=\frac{1}{3}\sum _{k}\left(\frac{{[{\rm{A}}]}_{j,k}}{{[{\rm{B}}]}_{j,k}}\right)1$$
(2)
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}.
$${\bar{c}}_{j,k}=\frac{1}{S}{\sum }_{i}\frac{{c}_{i,j,k}}{\frac{1}{4}{\sum }_{j}{c}_{i,j,k}}$$
(3)
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).
$$\overline{{[{\rm{A}}]}_{j}/{\bar{c}}_{j}}1=\frac{1}{3}\sum _{k}\left({[{\rm{A}}]}_{j,k}{\bar{c}}_{j,k}\right)1$$
(4)
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^{52}. This spatial separation allows the entire chamber to be frozen. The content of the chamber is then cut into four equalsized 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^{2} s^{−1} for amino acids^{53} and D_{i,AA} ≈ 1,400 µm^{2} s^{−1} for the 2aminoazoles and nucleotide components^{54}. 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
$$\rho ({\bf{u}}\cdot \nabla ){\bf{u}}=\nabla \cdot \left[p{\bf{I}}+\mu (\nabla {\bf{u}}+{(\nabla {\bf{u}})}^{{\rm{T}}})\frac{2}{3}\mu (\nabla \cdot {\bf{u}}){\bf{I}}\right]+\rho {\bf{g}}$$
(5)
and the continuity equation in the steadystate case
$$\nabla \cdot \left(\rho {\bf{u}}\right)=0$$
(6)
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 timedependent drift–diffusion equation
$$\frac{\delta {c}_{i,{\rm{num}}}}{\delta t}=\nabla \cdot \left[{D}_{i}\nabla {c}_{i,{\rm{num}}}\left({\bf{u}}{\boldsymbol{}}{S}_{T,i}\cdot {D}_{i}\nabla T\right){c}_{i,{\rm{num}}}\right]$$
(7)
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 custommade LabVIEW program incorporating Levenberg–Marquardt algorithms (see SingleNTD_TrapFitter_V16.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). Steadystate concentration profiles shown in Extended Data Fig. 2c,d and Supplementary Fig. 6 were obtained by calculating
$${c}_{i}(\,y)=\exp \left(\frac{\frac{{q}_{i}}{120}}{1+\frac{{q}_{i}^{2}}{\mathrm{10,080}}}{S}_{T,i}\Delta T\frac{y}{\alpha }\right),\,\text{with}\,q\equiv \frac{\Delta T\beta g\rho {\alpha }^{3}}{6\eta {D}_{i}}$$
(8)
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 steadystate 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^{−1} (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 selfmade LabVIEW program that fills a 3D array (see NetworkSimulation_v31.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.
Glycinedimerization experiments
For the experiments shown in Fig. 5, we closely followed the protocol for the dimerization of glycine^{11}.
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_{1}–k_{5} of the simplified reaction scheme^{11,12}:
$${\rm{Gly}}+{\rm{TMP}}\mathop{\to }\limits_{{k}_{1}}{\rm{GlyAct}}$$
(9)
$${\rm{GlyAct}}\mathop{\to }\limits_{{k}_{2}}{\rm{Gly}}$$
(10)
$${\rm{GlyAct}}+{\rm{Gly}}\mathop{\to }\limits_{{k}_{3}}{\rm{GlyGly}}$$
(11)
$${\rm{TMP}}\mathop{\to }\limits_{{k}_{4}}{\rm{waste}}$$
(12)
$${\rm{GlyGly}}\mathop{\to }\limits_{{k}_{5}}{\rm{Gly}}+{\rm{Gly}}$$
(13)
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_{1}–k_{5}, 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_{4} 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_{1}–k_{3}, k_{5} 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_V46.llb):
$$\frac{{\rm{d}}\left[{\rm{Gly}}\right]}{{\rm{d}}t}={k}_{1}\times \left[{\rm{TMP}}\right]\times \left[{\rm{Gly}}\right]{k}_{3}\times \left[{\rm{GlyAct}}\right]\times \left[{\rm{Gly}}\right]+{k}_{2}\times \left[{\rm{GlyAct}}\right]+{k}_{5}\times 2\times [{\rm{GlyGly}}]$$
(14)
$$\frac{{\rm{d}}\left[{\rm{GlyAct}}\right]}{{\rm{d}}t}={k}_{2}\times \left[{\rm{GlyAct}}\right]+{k}_{1}\times \left[{\rm{Gly}}\right]\times \left[{\rm{TMP}}\right]{k}_{3}\times \left[{\rm{GlyAct}}\right]\times \left[{\rm{Gly}}\right]$$
(15)
$$\frac{{\rm{d}}\left[{\rm{GlyGly}}\right]}{{\rm{d}}t}=+{k}_{3}\times \left[{\rm{GlyAct}}\right]\times \left[{\rm{Gly}}\right]{k}_{5}\times [{\rm{GlyGly}}]$$
(16)
$$\frac{{\rm{d}}\left[{\rm{TMP}}\right]}{{\rm{d}}t}={k}_{1}\times \left[{\rm{Gly}}\right]\times \left[{\rm{TMP}}\right]{k}_{4}\times [{\rm{TMP}}]$$
(17)
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_{4} 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 custommade LabVIEW program using a Levenberg–Marquardt algorithm (see FreeFitter_V46.llb). The rates obtained in this way are: k_{4} = 3.5 × 10^{−7} s^{−1} (±19%) and k_{1} = 1.0 × 10^{−3} Ms^{−1} (±41%), k_{2} = 1.1 × 10^{−4} s^{−1} (±10%), k_{3} = 9.7 × 10^{−5} Ms^{−1} (±21%) and k_{5} = 1.2 × 10^{−6} s^{−1} (±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^{−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_V46.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^{−3} K^{−1}. 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^{−1} is expected to be on the order of Volume_{chamber}/Q_{flow} ≈ 10^{1}–10^{2} 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).
[ad_2]
Source link