Machine learning enhanced droplet microfluidics Open Access

The formation of droplets in microfluidic devices has been the subject of significant interest in the past few decades^1–3 for many applications such as single-cell analysis,⁴ spheroids and organoids production,⁵ and digital polymerase chain reaction (PCR).⁶ The droplet formation mechanism in micrometre-scale devices is a complex phenomenon controlled by several parameters,^3,7–11 including surface tension, viscosity of the two non-miscible liquids, surface wetting properties, inertia, elasticity, addition of surfactant to either of the two liquid phases, and device geometry. Such complexity has so far prevented the establishment of a universal framework for the production of droplets in microfluidic devices, while the best approach is to empirically derive some scaling arguments that can only be applied within the experimental configurations presented in the study (see, for instance, the following articles^12,13 and references therein).

Machine learning (ML) and artificial intelligence (AI) have recently emerged as a very powerful tool to identify patterns across a complex set of data and to then deploy models to act on the data in a human-like manner but on a very short timescale.^14–16 ML and AI are sometimes used as synonyms, but they are different:^17,18 AI is capable of acting on a large set of data in order to trigger decisions on a short amount of time; ML is instead a class of AI-based methodologies that enable algorithms to learn without direct programing.¹⁵ For instance, Nitta et al.¹⁹ have employed ML to train the algorithm to recognize several populations of flowing cells in microfluidic devices; after recognition of the cell population type, the AI activated a deflector in the microfluidic device to sort cells according to a specific target.

Considering its power, it is not a surprise that ML has been used in the attempt of tackling the problem of droplet formation in microfluidic devices.^17,20 Most studies^21–25 have combined ML with the dimensional analysis in order to predict the experimental parameters required to produce droplets of a given size. Chagot et al.²⁶ employed ML to establish the effect of surfactant concentration on the size of water-in-oil droplets formed in flow-focusing configuration. ML has also been used in the context of droplet characterization,^27,28 sorting,²⁹ microrheology with optical tweezers,³⁰ bubble motion,³¹ droplet tracking,²⁷ computational fluid dynamics (CFD),³² and experimental optimization for digital microfluidic applications.³³ However, the droplet microfluidic studies introduced earlier have been limited in their analysis, most often not allowing the user to depart from spherical droplets nor enabling selection of the droplet generation rate, which is an essential parameter, especially in the context of encapsulation.³ Recently, Lashkaripour et al.³⁴ presented a detailed study where ML was employed to derive the dimensions of a flow-focusing microfluidic device and the experimental volumetric flow rate values required to produce droplets of a given size at a given generation rate, both pre-selected by the user. The authors also enabled the possibility to upload a new set of data to generate new predictions based on the requirements from the users, who may need to employ a combination of fluids different from the water and mineral oil presented in the study.

Despite the progresses made in the area of ML applied to droplet microfluidics, there are still some unaddressed problems. First, the majority of studies focused on the simple combination of water and mineral oil, which is not the preferred combination of liquids employed in biomedical applications, where fluorinated oil and phosphate buffer saline (PBS) are generally preferred.³⁵ Second, most of the studies relied on the possibility for the users to have access to microfluidic fabrication facilities, something that is not always true and that can hinder deployment of microfluidic solutions in remote areas. Finally, many of the papers published so far used and optimized a single architecture based on previous works or preference. A comparison between different approaches is essential, particularly for future work that may consider more parameters or more complex fluids.

In this work, we employed ML to predict the experimental parameters required for droplet generation in commercialized microfluidic devices using PBS and fluorinated oil as dispersed and continuous liquid phases, respectively. We compared the accuracy of three ML architectures in predicting the experimental conditions to obtain droplets of a given size and generation rate. We also provided a proof of concept related to the encapsulation of algae in both PBS and algae media. Finally, we developed a simple standalone app that can be used locally in conjunction with commercial microfluidic devices to predict the values of the volumetric flow rate of both the continuous and dispersed phases required to obtain droplets having a pre-determined size and generation rate.

We employed phosphate buffer saline (PBS, Sigma Aldrich UK) as the dispersed phase and fluorinated oil FluoSurf 2 wt. % in HFE 7500 (Dolomite Microfluidics, UK) as the continuous phase. We used a force tensiometer (Biolin scientific) to measure the interfacial tension γ between the FluoSurf and the PBS at room temperature (⁠ $T = 20$ °C), obtaining a value of $γ = 3.39 ± 1.09$ mN/m. The viscosity of the FluoSurf at $20$ °C was evaluated using a Twin Drive rheometer working in a strain-controlled mode (Anton Paar MCR702) equipped with a 50 mm cone and plate metal configuration, obtaining a constant value of $μ c = 1.56 mPa s$⁠.

We employed three commercially available microfluidic flow-focusing devices (Dolomite Microfluidics, UK) having three etching depths, namely, $100 μ$m (part 3200512), $190 μ$m (part 3200506), and $275 μ$m (part 3200824), all having a fluorophilic coating on their surfaces, meaning that they could be used for the production of aqueous droplets when the continuous phase was fluorophillic (Fig. 1).

The flow was controlled using a syringe pump having two independent flow units and also a pressure sensor on each flow stream (Dolomite Quad Pumps), equipped with the Mitos Quad Pump Green syringe pairs (Dolomite Microfluidics), with a lower flow rate limit of $5 μ$l/min and an upper limit of $1250 μ$ l/min. The advantage of this system was the possibility of generating a continuous flow of liquids taken from a reservoir at a constant volumetric flow rate (at variance with more traditional syringe pumps that require the syringes to be filled in as soon as the liquid in the syringe ends), while also simultaneously enabling the measurement of the overall pressure drop in the system. This approach was essential to monitor pressure drop variations and accumulation points of pressure that could bias the accuracy and reliability of the experimental data set.^36,37 The lower end of the volumetric flow rate was $5 μ$l/min for both the dispersed or continuous phases (this was set by the type of syringes employed), while the upper flow rate limit was reached when the overall pressure drop for either of the two flow streams reached 10 bars which, in our experiments, roughly corresponded to a volumetric flow rate of $120 μ$l/min. According to the manufacturer, the maximum pressure tolerable by the microfluidic device was 25 bars, meaning that when the volumetric flow rates of the two streams were equal to $120 μ$l/min, the overall pressure drop was around 20 bars, close to the device limit. In the case of overpressure, the flow was automatically arrested by the syringe pumps.

The connections between the microfluidic device and the syringe pumps were made of fluorinated ethylene propylene (FEP) commercially available tubing with an internal diameter of $100 μ$m and an external diameter of 1.6 mm (Dolomite Microfluidics). The choice of the internal diameter of the tube was driven by the requirement to obtain a stable flow at flow rate values of $5 μ$l/min, which was the lower end investigated here for both flow streams. We observed that the droplet formation at the junction for flow rate values of $5 μ$l/min was unstable when employing the same FEP tube having an internal diameter of either $800$ or $250 μ$m (data not showed). The reason for such a phenomenon was the fact that the overall pressure drop in the system (given by the sum of the pressure drops in the tubing and the device) was comparable with the pressure drop at the junction where droplet formation occurred, thus causing significant fluctuations in the flow, in line with well-established findings in the literature.³⁸ 

The flow in the microfluidic channel was observed using an inverted microscope (Zeiss Primovert) connected to a high-speed camera (Photron Mini Ux-50) with videos acquired at frame rate values varying between 125 and 8000 fps. The recorded videos were manually analyzed with droplets annotated using Matlab to derive the droplet length L, the droplet height H, and the generation frequency f; a minimum of 20 droplets per case was analyzed: this was deemed sufficient as the droplet size and frequency remained constant for the whole duration of the experiment, as also expected in droplet microfluidic applications.¹⁰ 

Three machine learning architectures, labeled as architecture 1, architecture 2, and architecture 3 were compared in this study (Fig. 2). Architecture 1 was made of a series of two neural networks with three hidden layers.

The layers consisted of a number of nodes and associated activation functions (in our case reLU) to determine the importance of a particular neuron when predicting the output. The layers were connected to each other via weights, which were subsequently updated during the learning process. The application of hidden layers enabled the network to account for dependency between variables along with complex relationships between input and output parameters.

Our overall challenge was summarized by a multi-output regression task, and therefore, the number of target variables was set at the number of nodes in the output layer. Hereafter, we will refer to the input parameters as those required by the user (e.g., droplet size), while the output parameters will always refer to the predictions made by the machine learning (e.g., volumetric flow rate values).

The operating principle of architecture 1 is described in Fig. 2(a). First, the droplet length L and the droplet height H were used as input parameters to then predict the volumetric flow rate values of both continuous (Q_c) and dispersed (Q_d) phases. After this first step, the values of Q_c and Q_d estimated from the first neural network, together with the droplet length L, were fed as inputs into a second neural network to predict the value of the droplet generation frequency f. In practical terms, architecture 1 enabled the prediction of Q_d, Q_c, and f, given an initial set of input related to droplet length L and droplet height H. This meant that the generation frequency f could not be specified by the user a priori, but it could only be the result of the neural network.

Architecture 2 took a sequential modeling approach based on standard vector regression (SVR), which practically means an architecture based on the fitting of the training data [Fig. 2(b)]. Via the use of kernels, the algorithm could fit a hyperplane in high-dimensional space [green line in Fig. 2(b)], which could be used to estimate the relationship between the input and output parameters. Architecture 2 was easily able to deal with correlations between output variables; however, its success was dependent on the order of chained regressions. Here, L and H were employed as input parameters (similarly to architecture 1), but only Q_d was predicted via the first SVR. Furthermore, the predicted value of Q_d was employed together with the original input parameters to predict Q_c using a second SVR. Finally, L, Q_c, and Q_d were used to predict the value of f. This architecture was very rudimentary and did not rely on neural networks.

Architecture 3 was a single neural network presenting L, H, and f as input parameters and Q_c and Q_d as output parameters [Fig. 2(c)]. Architecture 3 was the most flexible among the three with respect to droplet microfluidic applications, as it enabled the user to set a priori droplet length, height, and generation frequency. The reliability of each architecture with respect to either the ground truth data (droplet size) or the mass balance (generation frequency) was assessed using the coefficient of determination R².

The choice of these architectures rested on some considerations made in relation to the existing state of the art. Machine learning algorithms similar to architectures 1 and 3 have been used separately in different studies but, to the best of our knowledge, they have never been directly compared. Such comparison is important to elucidate whether training the algorithm with the droplet generation frequency is desirable or not. Additionally, we wanted to see whether a neural network was required at all to enhance droplet microfluidic application: for this reason, we tested architecture 2, which is based on an advanced fitting, rather than a neural network.

Tetraselmis suecica CCAP 66  $/$ 4 was obtained from the culture collection of algae and protozoa (CCAP) and cultivated in f $/$2 medium³⁹ prepared with seawater at $18$ °C and $80 μ mol m − 2 s − 1$ constant light and shaken at 90 rpm. We used the force tensiometer described earlier (Biolin scientific) to measure the interfacial tension γ between the FluoSurf and the algae medium at room temperature (⁠ $T = 20$ °C), obtaining a value of $γ = 8.69 ± 0.15$ mN/m.

The same experimental apparatus described before was used for the algae encapsulation experiments. For encapsulation in the algae media, the solution containing algae was employed as the disperse phase, and the experimental values of Q_c and Q_d were evaluated via architecture 3. For the experiments in phosphate buffer saline (PBS), the algae where first centrifuged using a micro-centrifuge, the media was removed using a pipette, and then the algae were re-suspended in the PBS, that was subsequently employed as the dispersed phase for the encapsulation experiments. The droplets generated were then collected on a Sedgewick rafter counting chamber, closed with a coverslip, and sealed within a petri-dish to prevent droplet evaporation.

To support this work, we designed a user-friendly graphical interface (GUI), enabling the user to easily predict the necessary flow conditions required to produce droplets of desired dimensions. Our app also enables the user to retrain our architectures using their own experimental data and use these to make predictions. The app was designed to tackle some of the common challenges associated with the task of droplet formation and, therefore, promote greater adoption within the microfluidic community. The GUI was developed using Rust, a multi-paradigm, general-purpose programing language that became prominent in the software industry in recent years due to the advantages of memory safety and better performance for certain tasks. All AI models were developed and trained using Tensorflow and Keras. The code for these steps was written in Python and called by our application using PyO3, a rust library that provides bindings for a Python interpreter.

The set of data employed to train a machine learning algorithm is generally called “ground truth” to make a clear distinction between the experimental or simulation data employed for the training, which are considered to be the “truth,” and the machine learning predictions. Therefore, it was critical that the ground truth was made of a highly reliable and reproducible dataset, or the predictions from the machine learning algorithm could be incorrect.^16,34 The formation of droplets in microfluidic devices is a very complex problem, which can be affected by several parameters, including pressure fluctuations, microfluidic fabrication, and surface wetting.^7–10 To minimize the impact of fabrication techniques and channel wetting on the droplet formation mechanism, we employed commercialized flow-focusing devices, whose reproducibility is ensured by industrial quality standards, while also reducing potential biases related to fabrication in different laboratory settings. Additionally, commercialized devices can be deployed in areas where microfluidic fabrication techniques, such as soft lithography,⁴⁰ milling,⁴¹ or 3D printing,^42,43 are not accessible.

Taken together, our results showed that the experimental dataset derived in the three microfluidic flow-focusing devices was robust enough to be employable as ground truth for the training of the machine learning algorithms. We were justified in employing a machine learning algorithm over a standard fitting procedure, as we could not clearly identify a mastercurve describing all our experimental results with the data following the scaling only qualitatively.

In this section, we analyzed the prediction efficacy of the three architectures with respect to the ground truth data and the theoretical scaling introduced in Sec. III A (Fig. 4).

First, we focused on the microfluidic device with $100 μ$m of etching depth and compared the flow rate ratio predicted by the machine learning algorithm against the ground truth data [Fig. 4(a)]. We employed N_GT = 79 experimental points to train the three architectures and observed that architecture 3 delivered the most accurate predictions, while architecture 2 provided the least accurate results. These results are not surprising considering that architecture 3 presented the droplet generation frequency f as input parameters, and f was linked to Q_d via Eq. (2). The addition of the experimental ground truth frequency in the model training facilitated a more accurate predictions of the flow rate ratio required to achieve droplets of a given size.

Architecture 2 performed poorly, and this was also expected as this approach was based on a chained regression. Errors made by the algorithm when predicting Q_d would essentially accumulate, leading to a further deviation than typically expected when predicting Q_c. The number of ground truth data N_GT = 79 was low for the variety of phenomena experimentally observed: specifically, we could clearly observe a transition from droplets being arranged on two lines [Fig. 5(a)] to droplets being arranged on a single line [Fig. 5(b)], meaning that the overall number of dataset available for each scenario was low, thus causing architecture 2 to return incorrect predictions.

We also compared the machine learning predictions with the ground truth data and with the scaling in Eq. (1) [Fig. 4(b)]. The purpose of this numerical test was to both assess how well each of the models had trained along with comparing the data derived with the scaling previously reported in the literature.¹² We observed that the predicted data followed closely the scaling and the ground truth, which was also expected considering the fact that the scaling was found to describe better the ground truth data for the $100 μ$m etched device. Furthermore, we compared the data of droplet generation frequency obtained by three architecture with the mass balance in Eq. (2) [Fig. 4(c)], observing a behavior similar to our previous experiments, with architecture 3 being more accurate than the others, and with architecture 2 performing the worst among the three (⁠ $R 2 < 0$ for architecture 2).

We repeated the same investigation for the $190$ [Figs. 4(d)–4(f)] and $275 μ$m [Figs. 4(g)–4(i)] etched devices. We again observed that architecture 3 was the best performing architecture. It is important to re-state here that the deviations observed between the empirical scaling and the machine learning predictions showed the need for a data driven approach, as the predictions as well varied from the theoretical scaling while being in close agreement with the ground truth data [Figs. 4(e) and 4(h)].

In summary, we observed that the predictions for architecture 3 were the most reliable with respect to the experimental ground truth data. This observation was rooted in the fact that architecture 3 employed the droplet generation frequency as an input parameter, thus somehow constraining the value of Q_d, and leading only to the prediction of Q_c. Architecture 2, based on chained regression, was found to perform the worst, and this was ascribed to the fact that this architecture was based on a chained approach with consequently larger errors. We were not too surprised by the supremacy of the neural network (architectures 1 and 3): due to the employment of a number of hidden layers, neural networks are often able to account for additional parameters that are not explicitly defined. Providing adequate steps are taken to prevent overfitting, the network is typically expected to perform well.

One of the key problems with supervised machine learning is overfitting,^47,48 where the experimental noise is computed within the architecture during the training, thus resulting in predictions that are only valid for specific experimental conditions and cannot be generalized to other systems, or even easily replicated across different laboratories because of experimental variability. This problem is particularly relevant in droplet microfluidics,¹⁰ where the resulting droplet size, shape, and frequency are dependent among a large variety of parameters.^7,8 We already took some mitigation steps to address this problem, especially by choosing commercial devices, which are generally produced at scale in a highly controlled manner in order to guarantee homogeneity. However, it is important to see how the machine learning architectures would respond to experimental variability, in order to gain confidence that our predictions are not the result of overfitting. To this aim, we carried out the following new experimental campaign, carried out on different days from the original ground truth dataset. Depending on the architecture, we supplied the input parameters and then derived the values of Q_d and Q_c required to generate the droplets; then, we applied the resulting values to our experiments and then compared the predictions with the actual data [Fig. 6 and Figs. S1 and S2)], thus simulating the real-life use of our algorithm.

We compared the predictions against the experimental data for the $190 μ$m etched device (Fig. 6). For each experiment (labeled with consecutive numbers on the x-axis), we compared the architecture value (white dashed bars) with the experimental observed value (red crossed bars), for the droplet length L [Fig. 6(a)], the droplet height H [Fig. 6(b)], and the droplet frequency f [Fig. 6(c)]. We did not compare the data with architecture 2, because it was the least performing, as already described in Sec. III B [Figs. 4(a)–4(c)]. For droplet length and height, we observed excellent agreement between predictions and observed data. For the frequency, instead, the minor disagreement observed for the frequency data was not uncommon in the literature, as other studies observed a similar trend.³⁴ Furthermore, it is well possible that not all the droplet sizes and frequency combinations were possible in our system, as we have constrained the microfluidic device geometry. It is, however, important to point out that the experimental parameters of Q_d and Q_c derived via architecture 3 were always returning value of the frequency f on the same order of magnitude compared to the observed values. We also repeated the same experimental campaign on the $100 μ$m etched device (Fig. S1) and on the $275 μ$m etched device (Fig. S2) and observed very good agreement. Note that we performed a larger number of experiments with architecture 3 to test the accuracy of the system in generating droplets having similar sizes but different generation frequencies.

Taken together, our results show that our training approach for three architecture does not suffer from overfitting, being the predictions reproducible on experiments performed on different days, thus accounting for some experimental variability. We also confirmed that architecture 3 is the best performing one for all the etched devices employed in this work.

The formation of droplets in microfluidic devices has widely been used in the context of encapsulation applications.^3,49,50 In this context, the frequency of droplet formation has become an essential parameter to know before performing the experiment, so that such frequency could be synchronized to the frequency of equally-spaced objects approaching the encapsulation area to overcome the Poisson stochastic encapsulation limit.^44,46,51,52 The droplet size, instead, is a crucial parameters for microfluidic applications related to encapsulation and co-encapsulation of cells⁵³ together with other cells or functionalized particles.⁵⁴ As a proof of concept featuring the encapsulation of biological objects, we here employed architecture 3 to predict the droplet size and frequency in the $100 μ$m etched device with the aim of encapsulating the Tetraselmis suecica algae (Fig. 7). We performed encapsulation experiments using both PBS [Fig. 7(a)] and Algae Media [Fig. 7(b)] as dispersed phase. The encapsulation of algae has been employed before to demonstrate intra-droplet growth,³⁵ as well as to elucidate inter-kingdom interactions in confined environments.⁵⁵ In our experiments, we clearly observed successful in-flow encapsulation at the flow-focusing microfluidic junction [Figs. 7(a) and 7(b)]. Notably, we compared architecture 3 values (the input parameters from which we obtained the values of Q_c and Q_d) and the ones observed in the experiments [Figs. 7(c)–7(e)]. For both PBS and Media, we observed good agreement on droplet sizes between the observed parameters and the input parameters in architecture 3. Most notably, we observed only a minor deviation when using the Algae media instead of PBS as the dispersed phase, despite the different value of interfacial tension measured for the two fluid combinations, in agreement with the literature.⁵⁶ Such observation can be justified in light of the weak dependency of the droplet length L with the Capillary number Ca, as observed from the ground truth data [Fig. 3(b)], with $L ∝ Ca − 0.35$⁠. For the frequency, while the observed value fell within a reasonable order of magnitude of the given values, there were more significant discrepancies, in line with similar type of deviations previously observed in the literature.³⁴ In other words, our results are expected to be usable also for other type of dispersed aqueous Newtonian phases which not depart significantly from the PBS, for instance with similar viscosity and interfacial tension values. We collected downstream the droplets with the encapsulated algae [Figs. 7(f) and 7(g)] and observed that the droplets retained their size and shape for both PBS and algae media.

The take-home point from this section was the fact that we could use architecture 3 for the purpose of algae encapsulation, even when the dispersed phase presented some different properties compared to the training fluid (PBS).

We developed an application based on architecture 3 using the Rust Programing Language (Fig. 8). The user can either use our build in model or insert their own trained model by activating the toggle on the top left of the panel: this is useful if the liquid combination departs significantly from the one employed in this work, especially in terms of viscosity and interfacial tension. On the top right, instead, the user is invited to define which etching depth the channel has, with choices among $100$⁠, $190$⁠, and $275 μ$m (which are the devices employed in this work). The users then insert the input data, namely, the droplet length L, the droplet height H, and the droplet generation frequency f. After inserting the input parameters, the app will return the predictions for the volumetric flow rate values of PBS (Q_d) and FluoSurf (Q_c). According to the results presented earlier in this manuscript, we expect the app to provide relatively accurate values even if the fluid combination is different from the one employed here, as long as the interfacial tensions between the two liquids is close to the 3 mN $·$ m (measured here), and the viscosity values for the two phases are close to the ones reported in our work.

In this work, we employed ML to predict the experimental parameters required for droplet generation in commercialized microfluidic devices using phosphate buffer saline (PBS) and fluorinated oil as dispersed and continuous liquid phases, respectively. We showed that the experimental dataset derived in the three microfluidic flow-focusing devices was robust enough to be employable as ground truth for the training of the machine learning algorithms. We were justified in employing a machine learning algorithm over a standard fitting procedure, as we could not clearly identify a mastercurve describing all our experimental results, with the data following the scaling only qualitatively. We observed that the predictions for architecture 3 were the most reliable with respect to the experimental ground truth data. This observation was rooted in the fact that architecture 3 employed droplet generation frequency as an input parameter, thus somehow constraining the value of Q_d, and leading only to the prediction of Q_c. Architecture 2 based on chained regression was found to perform the worst, and this was ascribed to the limited number of data available for the training. We performed a further experimental campaign demonstrate the validity of the predictions made by architecture 3 on different days, where experimental variation is a potential issue. We also demonstrated the use of architecture 3 predictions in the context of Algae encapsulation. Finally, we introduced an app where the users could select the device type, the droplet size, and generation rate to derive the predicted values of the liquid flow rates required.

We expect our results to be useful to the wider microfluidic community, being droplet microfluidic employed across a variety of different fields.^1–3 Future directions should include study of non-Newtonian droplet formation,^12,57,58 which has currently gained further attention in the context of controlled encapsulation.^44,46 Our results could also be employed to help tackling the problem of stochastic encapsulation, which prevent wide-spread adoption of such technologies within the medical setting.

See the supplementary material for Fig. S1 (predictions and experiments in the 100 μm channel) and Fig. S2 (predictions and experiments in the 275 μm channel).

F.D.G. acknowledges support from the EPSRC New Investigator Award (Grant No. EP/S036490/1) and from the Royal Society Research Grant (Grant No. RGS/R1/221263). F.D.G also acknowledges support via the Ser Cymru programme – Enhancing Competitiveness Equipment Awards 2022–2023 (Grant No. MA/VG/2715/22-PN47). This research contributes to the IMPACT operation, which has been part-funded by the European Regional Development Fund through the Welsh Government and Swansea University. We are grateful to Dolomite Microfluidics for providing the three microfluidic devices and the FluoSurf continuous phase free of charge.

The authors have no conflicts to disclose.

Claire Barnes: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – review & editing (supporting). Ashish Sonwane: Software (lead). Eva C. Sonnenschein: Data curation (supporting); Resources (supporting); Writing – review & editing (supporting). Francesco Del Giudice: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request. The application reported in Fig. 8 is available in GitHub at Ref. 59.