Development of a predictive response surface model for size of silver nanoparticles synthesized in a T-junction microfluidic device

Optimisation of the parameters governing the synthesis of silver nanoparticles (AgNPs) is a critical step in controlling particle size, which facilitates their application in diverse range of industrial and consumer related products. A T-junction microfluidic system was used together with design of experiments, regression-analysis and response surface methodology to build a predictive numerical model for the size of silver nanoparticles (AgNPs). Aqueous solutions of silver-precursor and reducing/stabilizing agent were supplied by two separate channels meeting at the T-junction, with the reaction occurring downstream in the outlet tube. To improve the mixing of the reagents, the output tube was coiled onto a 3D-printed helical shape device, exploiting the creation of Dean vortices. The effects of both reaction and hydrodynamic conditions including the solution pH, collection temperature, helical curvature, flow rates and concentration of stabilising agent were investigated using a D-optimal experimental design. The obtained experimental size distributions for the AgNPs were fitted to a polynomial model with an average prediction error of around 13% and a 37 % maximum error. The model predicted the optimal synthesis conditions for the global and local-minimum sizes of AgNPs with an error of around 7.0% and 16.1% respectively. The average prediction error of the testing set was estimated to be 6.8% with 16.1% being the maximum error.


Introduction
Silver nanoparticles (AgNPs) are very well-known for their unique antibacterial (Goel, 2021), antiviral (Saadh et al., 2021) and anticancer (Sukumar et al., 2022) activity, high electrical conductivity (Wang et al., 2021), optical (Thanh Nguyen et al., 2022;Wan et al., 2022;Yen et al., 2015) and catalytic properties (Huff et al., 2020;Sharma et al., 2022).These allow their application in a diverse range of industry sectors including food (Zhao et al., 2022), diagnostics (Hasanzadeh et al., 2019), pharmaceuticals (Gul, 2021), and electronics (Bariya et al., 2018).Particle size distribution and shape are key factors determining these properties which are largely dependent on the choice of fabrication method and the control of operating parameters (Nathanael et al., 2022).Despite the availability of robust chemical (Quintero-Quiroz et al., 2019;Mei et al., 2021) physical (El-Khatib et al., 2018;Menazea, 2020) and biological (Mathivanan et al., 2019) methodologies enabling size-controlled fabrication of AgNPs, there is an increasing demand for sustainable and reliable routes to decrease the environmental impact and transport limitations of current batch methods.Microfluidic based synthesis (Liu et al., 2018;Wojnicki et al., 2019) has gained significant interest as it permits chemical reactions to be carried out under highly controlled hydrodynamic conditions with reduced energy and material consumption in narrow channels of (10 − 4 -10 − 3 ) m.These types of reactor can be used in continuous flow (Zhu et al., 2021) and droplet-based (Xu et al., 2016) operation.However, continuous flow systems present some challenges for nanoparticles synthesis including poor radial mixing and limited axial mixing caused by laminar flow profiles (Gonidec and Puigmartí-Luis, 2019).These limitations can be overcome by the use of coiled and helical channels which promote radial mass transfer due to formation of Dean vortices whilst retaining narrow residence times (Gao et al., 2020).The formation of AgNPs involves three different steps, firstly the reduction of silver ions to silver atoms, secondly a nucleation step where the smallest thermodynamically stable clusters are formed and finally their growth to produce AgNPs (Thanh et al., 2014).
The design and optimisation of synthesis of AgNPs is an inherently complex process which involves many variables.For example, the degree of mixing of reagents (Kisyelova, et al., 2016), redox potentials of reducing agents (Demchenko et al., 2020), reaction temperature (Kaviya et al., 2011) and pH of the solution (Alqadi et al., 2014) are key parameters which can be tuned to control the reduction rate of silver ions and therefore the physicochemical characteristics of nanoparticles.Employment of one factor at time (OFAT) experimentation (Eckert and Trinh, 2016) to elucidate the effect of these variables individually whilst the others are held constant is a time and resource intensive process which ignores the effect of interactions between variables and results in the development of a poor prediction.To overcome these limitations, design of experiments (DoE) methods can be exploited (Mason et al., 2003).DoE methodology (Durakovic, 2017) is a statistical tool used to find the relationship between a set of independent variables (input variables) and a dependent variable (measured response variable) where factors are the input experimental parameters which can be quantitative or qualitative and levels are the different values that factors can take.
The most commonly used designs in scientific and engineering disciplines include full and fractional factorial (Altayb, 2021;Echeverría, 2021;Kania et al., 2021) and response surface designs (Nourafkan et al., 2017;Ayyubi et al., 2022;Ho et al., 2022).Full factorial design (FFD) (Antony, 2014b) contains all possible combinations of selected factors and levels.Although it provides the best insight into a system, FFD might generate a large set of experiments and be very expensive in terms of cost and time if the number of factors and levels are high.Conversely, fractional factorial design (Antony, 2014a) focusses on the main variables that affect the response (interactions are usually confounded) with the minimum possible number of runs (Mason et al., 2003).
Designs associated with response surface methodology (RSM) include central composite design (CCD), Box-Behnken design (BBD) and non-standard designs such as D-optimal design.These designs are used to provide an understanding of the system and the relationship between factors and responses and also to enable system optimization using minimum required experimental runs (Montgomery, 2013).CCD (Wagner et al., 2014) is based on factorial designs where the centre points are augmented with a group of axial points called star points to estimate the curvature, while in BBD (Das and Dewanjee, 2018), there are no axial points and the possible combinations of experimental points are placed at the centre points of the edges of experimental domain.Doptimal design (Montgomery, 2013) is an alternative to consider in cases where there is an irregular experimental region, a non-standard model that deviates from the usual first or second order ones or an unusual sample size requirement.
The application of DoE has been expanded to many industry sectors including pharmaceutical, (Ho et al., 2022;Baghaei, 2021;Baroutaji et al., 2019) and food (Encina-Zelada et al., 2019;Santos et al., 2018) as a part of design, development and optimisation of products and processes.There are a few studies to optimize AgNPs formulations using a systematic experimental design with main focus on the control of reactants used.Prema, (Prema, 2022) optimized batch synthesis of AgNPs using a central composite design (CCD) with three independent variables (AgN0 3 concentration, green tea extract concentration and temperature) at three different levels.The dependent variables were the size of AgNPs, zeta potential and polydispersity index.A second order polynomial quadratic equation was used to describe the relationship between the different factors across each response and used to find the optimised conditions (1 mM AgN0 3 , 0.5% green tea extract and 80 • C).Sarkar, (Sarkar et al., 2021) optimized the yield of AgNPs produced from Eucalyptus globulus fruit extract using a batch process.The authors used factorial design to identify the significant parameters and a CCD with three factors (AgN0 3 concentration, dose of the extract, time) at two levels.They found that the optimal conditions for maximum yield of AgNPs are AgN0 3 1.8 mM, 10% fruit extract and 12 h.
Whilst batch processes to produce AgNPs have been optimised using DoE approaches, continuous microfluidic approaches have not received the same attention despite their several benefits.Microreactors which are high-throughput experimental methods allow the use of less volume of reagents contributing in much more sustainable production to reduce the ecological footprint.DoE optimization is advantageous over other machine learning methods for microfluidic synthesis of AgNPs because it is an emerging approach where datasets are limited.This study investigates the use of an experimental design strategy based on the Doptimal approach to build a model to optimise the synthesis of AgNPs in a T-junction microfluidic device.The effect of chemical and thermodynamic parameters (solution pH, relative concentration of stabilising agent, collection, and storage temperature) and also parameters affecting the hydrodynamics (shape of reactor and reagent flow rates) on the size of AgNPs are investigated.This study aims to tackle the limitations of traditional OFAT, and batch-style chemical reactors approaches which are characterised by high energy and material consumption and impossibility of accounting for interactions between factors to deliver a more sustainable and controlled AgNPs synthesis route.Considering the increasing interest in nanomaterials and sustainable production, the proposed paper focusing on microfluidic synthesis and DoE, could be useful to reduce the ecological footprint and avoid large volumes of reagents, repeated and extensive experimental tests.

AgNPs synthesis and characterisation
AgNP synthesis was carried out by the reduction of silver nitrate in the presence of tannic acid, which has a weak reducing ability and trisodium citrate, which has a dual role in the reaction as a weak reducing and strong stabilizing agent.The concentrations of silver nitrate (SN) and tannic acid (TA) were 0.92 mM and 0.123 mM respectively, whereas concentration of trisodium citrate (TC) was varied in the range of 1.91-3.82mM.This protocol was based on previous work by Kašpar, (Kašpar et al., 2019) with modifications in the concentrations of TA and TC.Despite being smaller than the concentration of SN, the concentration of TA was chosen in excess, because each molecule of TA can reduce up to 20 molecules of SN (Sivaraman et al., 2009).The pH of the aqueous solution of tannic acid and trisodium citrate was set to either 7 or using 0.1 M NaOH.Silver nitrate, 99+% (AgNO 3 ), tannic acid (C 76 H 52 O 46 ), trisodium citrate dihydrate (Na 3 C 6 H 5 O 7 ⋅2H 2 O) and sodium hydroxide (NaOH) were purchased from Alfa Aesar.The microfluidic device was composed of 0.5 mm inner diameter polytetrafluoroethylene (PTFE) cylindrical tubing (Cole-Parmer).
The aqueous solutions of reagents were supplied by two separate inlet channels meeting at a Tee tubing junction (0.020 ′′ Thru-Hole) (Upchurch Scientific), whereas the reaction was carried out in the outlet channel of 2 m length as illustrated in Fig. 1.The inlet flow rates were controlled using syringe pumps (World precision instrument-AL-4000) equipped with 5 mL plastic syringes (Fisher) and were set to be identical for each inlet.The flow rate referred to in the below discussion is defined as the flow rate for one inlet.
The outlet tubing was either straight or coiled using 3D printed helix devices with a diameter of 3 mm.Particle size and zeta-potential was measured in triplicate from three independent experiments using Zetasizer Nano series, (Malvern Panalytical, Malvern, UK) with the refractive index of AgNPs set to 0.135.UV-vis spectroscopy measurements at nm which is characteristic for AgNPs (Shaikhaldein et al., 2020;Njagi et al., 2011) have been performed to indicate the presence of material obtained in the product stream.The results from dynamic light scattering (DLS) were confirmed using transmission electron microscopy (TEM) (JEOL JEM-1400).One drop of the sample was placed onto formvar/carbon 400 mesh copper grids (Agar scientific), followed by airdrying at room temperature for a few minutes.Particle size distribution was measured from the TEM images using Image J. The concentration of silver ions was measured using a silver electrode (EDT direction) connected to mV meter (Mettler Toledo-FP20) calibrated by standard solutions.Each of the samples was measured in triplicate.The final concentration of silver ions was tested for two different flow rates with two different outlet tube lengths to determine the reaction completeness.

Process optimization using DoE
The aim of the optimisation was to identify experimental conditions which produce the minimum size of nanoparticles.The preliminary experiments along with analysis of literature data from Marciniak, (Marciniak et al., 2020),Wu, (Wu et al., 2017); Das, (Das et al., 2021); Izak-Nau, (Izak-Nau et al., 2015) were used to determine the factors and levels for the experimental design.The optimal conditions for synthesis of AgNPs in microfluidic reactor were determined by response surface methodology and particularly, a D-optimal design (Montgomery, 2013) formed by five factors.Trisodium citrate concentration, flow rates, pH of the solution, curvature of outlet tubing, and storage temperature (numeric variables) were set at different levels, as shown in Table 1, while the size of AgNPs from DLS was investigated as the response.A Doptimal design was selected to limit the number of runs with a specific goal to minimise the variance of model parameters.It was preferred over other designs because there are multiple levels of factors, and the midpoint of levels is not exactly the central one.The final composite matrix was constituted of 162 experiments including three replications.The large number of experiments contributed to a higher accuracy in the predictions.
After the selection of representative runs for the experiment, regression-analysis was performed using Minitab and MATLAB.A quantitative dependence of the response on factors and interactions in the form of a regression equation was obtained using the least squares method.Analysis of variance (ANOVA) and lack-of-fit test were conducted to determine the statistical significance and quality of fit of the regression model at a confidence interval (CI) of 95% (Eq.S10).The goodness of fit for the model was expressed through the coefficient of determination (R 2 ), (Eq.S4) and the adjusted coefficient of determination (R 2 Adj), (Eq.S5).Both terms were used to measure the performance of the model while R 2 Adj provided more precision and reliability by considering the number of independent variables in the model (Cheng et al., 2014;Cho et al., 2019).Response surface plots were also used to visualise the relationship between the independent and dependent variables.A p-value ≤ 0.05 was considered statistically significant.Additional details may be found in the Supplementary Information.

Development of the model and statistical analysis
The experimental data (see Table S1 in SI) collected through the Doptimal design, were fitted to separate polynomial functions to predict the size of AgNPs (y) in nm, where Model 2: (2)  Model 3: (3) Different model structures were considered, however, examples of only three models are presented here by Eqs.1-3.Specifically, Eq. ( 1) shows all possible squared and cubic terms, Eq. ( 2) includes only all squared interactions, whilst Eq. (3) presents only some of the possible square and cubic terms such as The coefficient of determination (R 2 ), Eq.S4 and the adjusted determination coefficient (R 2 Adj), Eq.S5, for the models Eqs 1-3 (Table 2) are very close to 1 indicating that models provide a good fit to the data.However, R 2 and R 2 Adj are not reliable indicators for the formation of a good regression model.The value of R 2 increases as a new variable is added, even if it is not statistically significant.The R 2 Adj includes the number of parameters but it gives only the variation of independent variables which affect the dependent variable.The selection of the relevant factors in a regression model and consequently the number of model parameters, affects the complexity of the model and ultimately the Akaike Information Criterion (AIC) (Akaike, 1974) and Bayesian Information Criterion (BIC).Both AIC and BIC are widely applied in the area of machine learning and data science to discriminate between the accuracy of models by the number of parameters involved and thus identify those that best describe a system.In BIC, the penalty term is larger than that in AIC.Therefore, BIC generally results in choosing models with less parameters while AIC has larger chances of choosing a model with too many parameters (Claeskens and Hjort, 2008).The expressions of the statistical approaches used to select model are as follows (Akaike, 1974;Schwarz, 1978;Bhattacharya et al., 2016): where n is the number of experiments, k the number of parameters (or coefficients) in the regression model and SSE is the sum of squared errors, (Eq.S3).The model with smallest AIC and BIC is thus the optimal one.On the basis of the above criteria, Eq. (3) was selected over the others because it had the smallest values for both AIC and BIC (Table 2).This model can predict the size of AgNPs within the limits of the experimental factors with a minimum standard error equal to 1.84 nm.Note, the model (Eq.1) includes all variables and their interactions and therefore has the highest values of R 2 and R 2 Adj, but has larger values of AIC and BIC than model (Eq.(3) because of larger k.
The normal probability plot is a tool to determine if a data set follows a hypothesized distribution (vertical axis is adapted for normal distribution).The probability plot in Fig. 2A revealed that the plotted points lie on a straight line while histogram plot in Fig. 2B showed that variation is normally distributed.Both residual plots indicated that there is no deviation from normality, so the errors appear to be independent.
ANOVA (King et al., 2010) and luck-of-fit test (Fisher, 1922) were performed on the regression model (Eq.3) for the size of AgNPs to determine the statistical significance and the quality of fit of the model to the dataset.ANOVA uses the sum of squares Eqs.S1 -S3 to compare response and error data, the mean square Eqs.S6, S7 to provide an estimate of population variance and the F-test (Eq.S8) to show the significant effects at a given probability level based on the sample variance (Wagner et al., 2014).The high F-values along with p-values smaller than 0.05 obtained in the ANOVA test for the model showed the significance of the correlations between the experimental set of data and the fitted polynomial model.The results of the ANOVA and lack-of-fit tests are reported in detail in Table 3B.As the assumptions on the normality and variance of residuals were met, confidence intervals (CI) (Eq.S10) for regression parameters (Table 3A) were calculated.These upper and lower bounds for the CI showed the uncertainty for the parameters (or coefficients) with a confidence of 95%.The t-values (Eq.S9) which are also shown in Table 3A are used to establish the uncertainty associated with the estimates.A critical t-value (tcrit) equal to 1.98 was used to show how well each parameter was estimated (t > tcrit means that the factor is statistically significant) assuming the parameters are normally distributed.
The Pareto chart in Fig. 3 shows the standardized effects of factors, from the largest to the smallest, with a reference line to note which are statistically significant (Ramachandran et al., 2021).This chart was built based on t-values (Eq.S9) with the threshold calculated using the Student's t-distribution (Student, 1908) with N-Ntheta degrees of freedom (N is the overall number of points used in calibration, Ntheta is the number of parameters or coefficients).Using a pareto chart, the importance and magnitude of the factors on the size of AgNPs can be determined.Particularly, the most significant factors of the model (Eq. (3) include the pH of the solution, the concentration of trisodium citrate, the interaction of pH and trisodium citrate and the curvature of the output channel.All the experimental factors cross the reference line that is at 1.98 which means all of them are statistically significant at the 0.05 level with current model terms.

Response surface methodology
A graphical representation of the regression model obtained (Eq.( 3) is shown via the two-dimensional contour plots in Fig. 4. The plots were used to study the combined effect between the formulation parameters and the size of AgNPs and to define the optimum condition of each factor to deliver a minimum size of AgNPs.
The plot of storage temperature versus trisodium citrate concentration (Fig. 4B), shows that the size of the AgNPs is virtually independent of storage temperature for a trisodium citrate concentration below 3 mM.The AgNPs produced are also small under these conditions.There is an effect of storage temperature and concentration of TC on the size of particles above TC = 3 mM.The plot of curvature of the outlet channel versus storage temperature (Fig. 4D) shows that there is a strong effect on the size of particles at higher storage temperatures and smaller curvature of the outlet channel.According to the plot of storage temperature versus pH in Fig. 4E, the interaction of pH and storage temperature strongly influences the size of particles only at smaller pH.Similarly, the plot of curvature of the channel versus pH (Fig. 4C) of the solution shows a significant effect of curvature of the channel at a value of pH below 11.
The effect of different factors on the size of AgNPs, in particular Fig. 4A, will be discussed in more detail in the next sub-sections.

Effect of trisodium citrate concentration, pH and storage temperature
Trisodium citrate acts as both reducing and stabilizing agent in the synthesis of AgNPs.Citrate ions can reduce silver ions at elevated temperatures (Khatoon et al., 2017) or under alkaline conditions (Marciniak et al., 2020), but they can also weakly interact with silver generating a charged layer that works as an electrostatic barrier to aggregation (El Badawy et al., 2012).The effect of trisodium citrate concentration was investigated over a range from 1.91 mM to 3.82 mM.Fig. 4A shows that there is an optimal range of concentration of TC for the minimum size of AgNPs between 2.3 mM and 3.1 mM at pH 7. It may be assumed that trisodium citrate within this concentration range fully covers the surface of the AgNPs and thus prevents further particle growth.At lower concentrations, the surface coverage is incomplete enabling the particle growth.An increase of concentration above 3.1 mM can result in a decrease of Debye length due to an increase of ionic strength of the solution, below a certain critical point aggregation may occur.Henglein and Giersig (Henglein and Giersig, 1999) investigated the formation of AgNPs through γ-irradiation of AgClO 4 solution which contains 2-propanol and a range of citrate concentrations from 0.05 to 1.5 mM.In this system, citrate did not participate in silver reduction and so is present solely to stabilise the AgNPs.It was found in (Henglein and Giersig, 1999) that the small particles with a narrow size distribution were found for a range of citrate concentrations between 0.1 and 0.5 mM.At smaller concentrations, large particles with a broad size distribution were formed.This was ascribed to a rather loose citrate layer on the surface of Ag clusters which was unable to prevent their growth first via condensation and later via coalescence.At concentrations larger than the optimum, larger particles were formed via coalescence due to high ionic strength of the solution.Das, (Das et al., 2021) also showed that there is an optimal concentration of trisodium citrate in AgNPs synthesis.
It is well-known that the pH of the solution often influences the chemistry of a reaction mixture.Here, the effect of pH on the reaction was examined by adjusting the pH at two levels of 7 and 12.The nanoparticles formed at pH 12 and a concentration of TC of 1.91 mM were smaller compared to those produced at pH 7 and a concentration of TC of 1.91 mM, indicating that a basic pH altered the kinetics of reaction.(Fig. 4A).Under alkaline conditions, the dissociation of citrate and dicarboxyacetone (DCA), which is produced as a by-product of the decarboxylation of citrate during silver reduction, is shown in Fig. 5.The completely hydrolysed DCA 2− (Marciniak et al., 2020) and Cit 3− species (Dong et al., 2009) at pH 12 have a higher reducing ability, leading to faster reduction of the silver ions to silver atoms and thus smaller particles.
Furthermore, at basic pH, tannic acid can be hydrolysed into gallic acid and glucose (Osawa and Walsh, 1993) where the former is a good reducing agent and the latter is a good stabilizer.Hence, the chemical rearrangement in the structure of tannic acid when exposed to alkaline pH promotes the formation of smaller particles (La Spina et al., 2020).
However, when the pH of the solution was increased alongside a trisodium citrate concentration above 2.87 mM, the size of particles   increased (Fig. 4A).The balance between attractive van der Waals and repulsive electrostatic forces which governed the stability of colloidal dispersion can be manipulated by changes in pH, ionic strength and temperature (Chhabra and Basavaraj, 2019).In general, as the ionic strength of a solution increases, the thickness of double layer associated with the AgNPs is reduced, resulting in weaker electrostatic repulsions between the AgNPs (Phan and Haes, 2019).This can be the reason for an increase in particle size at high concentration of trisodium citrate.The surface charge of AgNPs is governed by the ionization of carboxyl-groups, ( − COOH), of the citrate molecules, which is dependent on the pH and ionic strength of surrounding aqueous phase (Badawy et al., 2010).AgNPs obtained at concentrations of TC less than 2.87 mM were smaller for both values of pH studied, when compared to AgNPs prepared at larger concentrations (3.82 mM) (Fig. 4A).At higher concentrations of TC, particles of smaller size were observed only at smaller pH.In Fig. 6A, it can be seen that the zeta-potential of AgNPs prepared at pH 7 and TC = 2.87 mM is approximately the same as the AgNPs obtained at pH 12 and TC = 1.91 mM, indicating that higher concentration of stabilizing agent increases the stability of AgNPs in the same way as alkaline pH.By increasing the pH, the carboxyl-groups ( − COOH) on the AgNP surface reach a higher degree of deprotonation and the value of negative surface charge increases.Thus, the electrostatic repulsion between particles is improved causing smaller AgNP diameters.On the other hand, AgNPs prepared at 3.82 mM TC showed a small gradual increase in diameter and a decrease in zeta potential absolute values with an increase in pH (Fig. 6A), indicating reduction in electrostatic repulsion due to decrease of surface potential and probably due to decrease of Debye length.the effect of pH and ionic strength on the growth of tannic acid-poly(Nvinylpyrrolidone) (TA-PVPON) multilayers, and found that the thicker multilayers are obtained at alkaline pH in the presence of high salt concentrations.They claimed that these conditions cause thicker multilayers due to screening of the negatively charged tannic acid.
The ageing of AgNPs dispersions can be also affected by storage temperature, since chemical reaction rate increases with temperature.
As it can be seen in Fig. 4E and D, the collection of samples at T = 0 • C quenched the reduction reaction or aggregation and led to the production of smaller particles.Many previous studies reported changes in nanoparticles samples kept at room temperature such as Velgosova, (Velgosova et al., 2017) and Peng, (Peng et al., 2010).

Effect of curvature and flow rates
In microreactors, variation in the flow rates of reagents and the flow configurations changes the residence times and mixing conditions, thus influencing the overall reaction of AgNPs including nucleation and growth steps (Baber et al., 2017;Ravi Kumar et al., 2012).The formation of AgNPs requires uniform reaction conditions because the reduction of silver ions to silver atoms is very fast.It has been reported that reduction stage can be carried out in less than 200 ms (Polte et al., 2012).In the present study, the flow rates of reagents varied from 1×10 − 6 dm 3 /min to 4×10 − 5 dm 3 /min and mixing was intensified by coiling the outlet tube on a 3D-printed helical design.The Dean number (De) (Eq.6), which represents the Dean force due to secondary flow in curved channels (Nivedita et al., 2017) was calculated to be 0.033, 0.34 and 1.39 respectively for the three different flow rates used in the experiments.
In Eq. ( 6), d is the inner diameter of the microchannel, R c is the radius of helix curvature and Re is the Reynolds number where ρ is the fluid density, u is the velocity of the fluid, and μ is the viscosity of the fluid.The increase in the flow rate decreases the residence time in the outlet tube, but increases the mixing, especially in helical channel.To check for completeness of reaction, experiments (pH 7, TC 1.91 mM, ST 0 • C, CR 0 dm − 1 ) were performed at two different flow rates (2 × 10 − 4 dm 3 /min and 6.7 × 10 − 5 dm 3 /min) with different outlet tube lengths (2 m and 0.667 m).The concentration of silver nitrate was measured at the exit using the silver electrode.Results showed that mixing (even in straight output channel) and residence time are both significant.At conditions of 2 × 10 − 4 dm 3 /min and 2 m length, the residual concentration of silver nitrate was 0.031 ± 0.0036 mM ( 3% of initial concentration) but at the same residence time at 6.7 × 10 − 5 dm 3 /min and 0.667 m length the silver nitrate concentration was 0.051 ± 0.0054 mM ( 5% of initial concentration).The larger residual concentration of silver nitrate in the last case shows that at the smaller flow rate, the reaction was farther from completion.Considering that reagent transfer in a straight tube is mostly diffusional, the degree of mixing at the T-junction, which increases with an increase of flow rate, may be a significant factor.The characteristic diffusion time of AgNO 3 was calculated as 36 s while the residence time was estimated as 117.75 s.At 6.7 × 10 − 5 dm 3 / min with a 2 m outlet tube, the residual concentration of silver nitrate was around 0.017 ± 0.00053 mM (less than 2% of initial concentration).Therefore, the reduction of silver nitrate was practically complete when the residence time was increased from 117.75 sec to 353 sec.It was also observed that the size of AgNPs is greatly influenced by the curvature of the output channel (e.g Fig. 4C).In straight channels under laminar flow, there is a characteristic parabolic profile.The curvature of the outlet tube modified the flow pattern by the formation of secondary flows (Dean vortices) which improved the radial mixing in the channel (Wu et al., 2017).Sufficient mixing allowed fast mass transfer during reduction process and homogenous growth and hence led to the formation of smaller particles.

Experimental validation of the model
The validation of the obtained model and verification of the optimization results was carried out according to the optimal synthesis conditions.Using (Eq.3), the model predicted the global minimum size to be 2.62 nm and a local minimum size of AgNPs to be 4.32 nm when the optimal conditions were set to pH 7, trisodium citrate 2.68 mM, curvature 33.3 dm − 1 , storage temperature at 0 • C and flow rate 4×10 − 5 dm 3 /min and to pH 12, trisodium citrate concentration 1.91 mM, curvature 33 dm − 1 , flow rate 2×10 − 5 dm 3 /min and storage of samples at 0 • C respectively.These predicted conditions for minimum size of AgNPs were confirmed experimentally with an error of 16.1% and 7.0% respectively, see Fig. 7.The size of AgNPs was measured for different conditions using dynamic light scattering and transmission electron microscope.A good agreement was found between the two characterisation methods (see Fig. S1 in SI).
More subsequent experiments were performed in triplicate to validate the predicted model.The predicted versus experimental values of AgNPs size for the training data and testing data are presented in Fig. 8A and B respectively.The mean absolute percentage error (MAPE) (Dash et al., 2017) was calculated using the equation below where y k is the measured response and ŷk is the predicted value: The training data contains 162 experiments, including replications, while testing data comprised of only 45 experiments including replications.The testing data were taken from independent experiments selected randomly within the limits of experimental factors.These include experiments close to optimal conditions and beyond them.The optimal points are illustrated with red dots in Fig. 8B.The confidence intervals (CI), Eq.S10, and prediction intervals (PI), Eq.S11, for the predicted responses were estimated to present the uncertainty with a probability of 95%.To be more precise, PI appears wider than CI because it shows the variability related to single observations, while CI reflects the variability associated with the uncertain estimates of parameters (or coefficients) (Lane and DuMouchel, 1994).
The predicted values are close to experimental responses of the testing data set and within the confidence and prediction intervals, suggesting that the model is successfully developed with correlation between the different factors and the size of AgNPs.The average prediction error of the training set was estimated around 13%, with 37% the maximum error while the average prediction error of the testing set was calculated to be 6.8% with 16.1% being the maximum error.

Conclusions
A D-optimal experimental design was applied to build a predictive regression model using response surface methodology to optimize the formation of AgNPs in a microfluidic continuous flow reactor.The statistical analysis provided an insight into factors effects on the measured response as well as interaction effects between factors which are difficult to identify using OFAT methods.Different response surface models were considered, and model discrimination was carried out to select the model using AIC and BIC.Through the investigation of five experimental factors including storage temperature, pH of the solution, concentration of stabilising agent (trisodium citrate), curvature of output channel and flow rates, it was shown that all factors are significant for the measured response while pH, concentration of trisodium citrate, curvature, and the interaction of trisodium citrate and pH are the most vital.The resulting model equation was graphically represented as a response surface function and allowed identification of optimized hydrodynamic and reaction conditions.The optimal combinations (pH 7, TC 2.68 mM, CR 33.3 dm − 1 , ST 0 • C and FR 4×10 − 5 dm 3 /min) and (2 × 10 − 5 dm 3 / min, pH 12, TC 1.91 mM, CR 33.3 dm − 1 and ST 0 • C) predicted AgNPs with an average size of 2.62 nm and 4.32 nm respectively.The average prediction error of testing data set was calculated as 6.8% confirming that a reliable model has been developed which correlates the effect of the different factors on the size of AgNPs.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Residual plots for the size of AgNPs A) probability and B) histogram of residuals.The data set follows a basic normal (Gaussian) distribution.
Fig. 6B also shows that as the concentration of trisodium citrate increases at pH 12, the surface-potential decreases.In line with this suggestion, Liu, (Liu et al., 2014) examined

Fig. 5 .
Fig. 5. A) Dissociation species of citrate, B) Production of dicarboxyacetone (DCA) as a by-product of the decarboxylation of citrate during silver reduction and C) Dissociation species of DCA.Reproduced from Gao and Torrente-Murciano (Gao and Torrente-Murciano, 2020) and Marciniak, Nowak (Marciniak et al., 2020).

Fig. 6 .
Fig. 6.A) Dependence of zeta potential of AgNPs on pH at different concentrations of trisodium citrate and B) Zeta potential of AgNPs at pH 12 for different concentrations of trisodium citrate.

Fig. 8 .
Fig. 8. Predicted and true particle size of the A) training data (162 experiments including replications) and B) testing data (45 experiments including replications) using Eq.(3).The confidence intervals (CI) and prediction intervals (PI) for the predicted responses were estimated with error 95%.

Table 1
Factors and levels for D-optimal design.
K. Nathanael et al.

Table 2
Values of R 2 , R 2 Adj, AIC and BIC for each model.

Table 3 A
) Parameters and B) analysis of variance for the fitted polynomial model (Eq.3) of AgNPs optimization.