Flexibility analysis for continuous ibuprofen manufacturing processes

Wenhui Yang ,Haoyu Yin ,Zhihong Yuan,*,Bingzhen Chen

1 State Key Laboratory of Chemical Engineering,Department of Chemical Engineering,Tsinghua University,Beijing 100084,China

2 School of Material Science and Engineering,Tsinghua University,Beijing 100084,China

3 Department of Chemical Engineering,Tsinghua University,Beijing 100084,China

Keywords:Continuous manufacturing Pharmaceuticals Mathematical modeling Flexibility analysis Mixed integer programing Optimization

ABSTRACT Continuous ibuprofen (a widespread used analgesic drug) manufacturing is full of superiorities and is a fertile field both in industry and academia since it can not only effectively treat rheumatic and other chronic and painful diseases,but also shows great potential in dental diseases.As one of central elements of operability analysis,flexibility analysis is in charge of the quantitative assessment of the capability to guarantee the feasible operation in face of variations on uncertain parameters.In this paper,we focus on the flexibility index calculation for the continuous ibuprofen manufacturing process.We update existing state-of-the-art formulations,which traditionally lead to the max-max-max optimization problem,to approach the calculation of the flexibility index with a favorable manner.Advantages regarding the size of the mathematical model and the computational CPU time of the modified method are examined by four cases.In addition to identifying the flexibility index without any consideration of control variables,we also investigate the effects of different combinations of control variables on the flexibility property to reveal the benefits from taking recourse actions into account.Results from systematic investigations are expected to provide a solid basis for the further control system design and optimal operation of continuous ibuprofen manufacturing.

1.Introduction

Ibuprofen is an analgesic drug with better efficacy compared to aspirin.Early clinical trials showed that 600 mg·d-1ibuprofen had the same curative effect as that of 3.6 g·d-1aspirin or 3 g·d-1paracetamol and high dosages up to 1000 mg·d-1showed no significant side-effect [1].Clinical studies of ibuprofen proved that it was not only effective for treating rheumatic and other chronic and painful diseases,but also showed great potential in dental diseases[2].Recently a report estimated that over 54 millions people use ibuprofen and its annual sales amount reaches 3 billion USD worldwide [3].In conclusion,ibuprofen is indeed a very useful compound.

As illustrated by Fig.1,Boots/Hoechst-Celanese(BHC)process is the epitome for synthesizing anti-inflammatory drugs like ibuprofen with batch production mode [4,5].Note that although batch manufacturing is a mature and dominant technology in pharmaceutical industry,it remains many endogenous problems.First,understanding of batch processes is still poor,thus operations are time-dependent and scale-dependent.Second,batch process with more stages means less yield,for example,a 6-stage Boots process has a lower overall conversion rate than that of 3-stage BHC process[4].Third,the storage of intermediates produced in each batch step costs a lot.Fourth,transfer efficiency of heat,mass and momentum limits the scaling up capacity of batch process [6].

In contrast with the batch manufacturing,continuous flow technology can provide a productive,green,and safe approach for ibuprofen manufacturing[7].Advantages of continuous synthesis over batch mode can be briefly expressed as follows.First,compared with batch process,continuous process controls temperature more accurately and transfers mass/heat more efficiently,which makes great contributions to the progress in yield [6,8].Second,continuous flow technology increases production capacity by adding the number of continuous lines in parallel instead of enlarging the scale of reactors so that the basic scale remains unchanged[9].Numbering-up instead of scaling-up avoids the limitation of transfer efficiency at the source and insures safety as well[10,11].Third,continuous manufacturing can reduce the waste,energy,and cost because of higher yield under better controlling and transferring condition.Clearly,continuous flow technology shows excellent potential in ibuprofen manufacturing.The patent from SRI International (formerly Stanford Research Institute) for the continuous flow driven ibuprofen manufacturing shows that the overall conversion rate can be up to 82% without much waste (Boots process has) or heavy and expensive catalyst (BHC process does) [4].

Fig.1.The BHC process for the batch ibuprofen manufacturing.

To be honest,up to now continuous ibuprofen manufacturing still stays in lab.Definitely,good process design and optimization from process systems engineering (PSE) perspectives can accelerate the development and deployment of continuous ibuprofen manufacturing technology.For instance,the simultaneous considerations of process design and process operability analysis can not only ensure the desirable dynamic operation performance,but also guarantee the drug quality.Operability analysis traditionally contains flexibility,controllability,reliability,and safety [12].For the ibuprofen manufacturing,flexibility analysis is in charge of the identification of feasible operation in face of variations on uncertain parameters.With regard to the continuous ibuprofen manufacturing process,uncertainties caused by mismatch on kinetic parameters and exogenous conditions are inevitable,it is therefore indispensable to execute the flexibility analysis to ensure feasible operations.The flexible region obtained by flexibility analysis is also referred to as the design space.It should be emphasized that flexibility analysis lays a solid foundation for other aspects of operability analysis.In fact,flexibility analysis has already been applied to the definition of design space in pharmaceutical industry,which has been proved very effective [13,14].

During the past several decades,flexibility analysis has been developed into one of relative mature branches to tackle uncertainties in the PSE community,and has already been applied into a great number of chemical processes such as heat exchanger networks synthesis[15,16],structural flexibility analysis for heat integrated distillation columns [17],and optimally designed flexible distillation sequences [18].Flexibility analysis was also launched on continuous processing plants -the multi-period optimization of practical refinery plant [19].Other applications include nature gas plant [20],multi-stage flash (MSF) desalination plant [21],and ternary distillation [22].All of these works agree well with a conclusion that with the help of flexibility analysis,the processes are able to adapt to the requirements of more uncertainties.

To further extend the scope of flexibility analysis,novel concepts and formulations have been proposed during the past few years.Ellipsoidal uncertainty set instead of traditional rectangle has been introduced to the flexibility analysis.For instance,Pulsipher and Zavala presented a mixed-integer formulation to compute the flexibility index under Gaussian uncertainty [23,24].Surrogate model-based approach for feasibility and flexibility analysis was proposed for the continuous tablets’ compaction process,which is essentially different from traditional mixed integer linear/nonlinear program (MILP/MINLP) based ones [25,26].The classification of uncertain parameters such as measured/unmeasured and exogenous/endogenous has also been paid with sufficient attention in recent papers [27,28].All of these excellent works,needless to say,excite the further theoretical and methodology development for flexibility analysis as well as the corresponding depth applications in the field of disruptive pharmaceutical manufacturing.

As aforementioned,the core of this paper is to investigate the flexibility property of continuous ibuprofen manufacturing while considering the computational issues associated with the flexibility index calculation.The major scientific contributions can be summarized as follows:

●We identify the flexibility index for the continuous ibuprofen manufacturing process with the consideration of uncertain process and model parameters.

●We provide a new reformulation for the flexibility index calculation and testify its validity and efficiency with several classical cases.

●We reveal the influences of different combinations of control variables on the flexibility index for the continuous ibuprofen manufacturing process.

The remainder of this paper is structured as follows: in Section 2,mathematical formulations of traditional flexibility analysis,and modified reformulations for flexibility analysis are presented.In Section 3,flexibility analysis is applied to continuous ibuprofen manufacturing process.Main concluded remarks are provided in Section 4.

2.Mathematical Formulations for Flexibility Analysis

2.1.Retrospective on the mathematical formulations

The optimization model for the flexibility analysis involves following variables and parameters.dis design variables,xis state variables,zis control variables,θ is uncertain parameters,θNis the nominal point inside the feasible region,Δθ-and Δθ+are the expected deviation of uncertain parameters,Tis the specified parameter ranges,δ is flexibility index,uis auxiliary variables.The involved equality and inequality constraints can be generically described as follows.

Clearly,equality constraints allow the elimination of state variables [12].

In a broad sense,the flexibility analysis problem can be divided into three categories: feasibility test problem,flexibility index problem,and optimal design problem [29,30].In this paper,to quantitively assess the flexibility property of the system,we mainly discuss the flexibility index problem.The flexibility index proposed by Swaney and Grossmann can be formulated as follows[12].

Note that there may be a confusing step here.In some cases,the formulations of flexibility index are maximum optimization problems.While they are minimum in other conditions,shown as Eq.(6).

It should be pointed out that the formulations represented by Eq.(5) and Eq.(6) are the same in essence.Fig.2 distinguishes the difference between them in the geometric sense.Eq.(5)intends to find the maximal hyper-volume inside the feasible region,while Eq.(6) intends to find the minimal hyper-volume on the boundary.In Fig.2,the dotted boundaries indicate the searching area and the purple boundaries are the identified optimal solutions.

Once introducing the auxiliary variableuand the Karush-Kuhn-Tucker(KKT)condition,Eq.(6)can be transformed to the following computationally tractable formulation Eq.(7).According to the theory developed by Grossmann and Floudas,the main sufficient conditions for solving Eq.(7) to global optimality are that constraint functions must be jointly quasi-concave inzand θ,and strictly quasi-convex inzfor fixed θ [31,32].

In Eq.(7),nz+1 stands for the number of active constraints,wherenzrepresents the dimension of the vector of control variablesz[32].Note that Eq.(7)takes control variables into consideration.Additionally,flexibility index calculation without control variables is represented by Eq.(8) [27].

2.2.Modified reformulations for flexibility analysis with uncertain process and model parameters

Measured parameters and unmeasured parameters mentioned in existing contributions [33,34] can be essentially referred to as process parameters and model parameters,respectively.Process parameters denote those directly measured ones,such as temperatures,pressures,concentrations,and flow rates.While model parameters are those priori unknown quantities in mathematical representation,and their values can only be estimated from experimental data sets.For examples,kinetic parameters,and heat transfer coefficients.The mathematical model for calculating flexibility index with process and model parameters while without any control variable can be presented as below [13].

According to Ref.[13],two available reformulations of model(9) exist: nested reformulation and decoupled reformulation.The logical structures of these two existing formulations are illustrated by Fig.3.The objective of both nested and decoupled methods is to reformulate the max-max-max trilevel problem into a single level problem with the help of KKT conditions.One level of the multilevel problem can be simplified as the solid and filled arrowhead shows.Detailed models for these two reformulations can be found in Supplementary Material.In this subsection,a new simplified reformulation will be proposed.

Fig.2.(a) The geometric sense of Eq.(5) and (b) the geometric sense of Eq.(6).

Fig.3.The logical structures of (a) nested reformulation and (b) decoupled reformulation.

Before providing our modifications on existing reformulations for calculating flexibility index,following lemma for those optimization problems without any inequality and equality constraints will be first proved.

Given.If the maximum off（x，y）exists,xandyare independent continuous variables,the following equality is true.

Proof.Ifxandyare two independent continuous variables.A necessary and sufficient condition expressed by Eq.(11)should be satisfied.

Then,an auxiliary variableuis introduced to represent the maximum value of the inner problem

KKT condition is then applied here to reduce the bilevel optimum problem into single level.The structure of the modified reformulation method is shown as Eq.(18) which can be used to calculate the flexibility index.

Therefore,the flexibility index problem can be formulated as the following MINLP formulation by the modified reformulation method.

The core essence of the simplified reformulation method is to merge the appropriate maximal optimal levels together.It should be emphasized that there is a significant assumption before using this modified reformulation method -the variables of θpand θmshould be independent and continuous.The first two level variables θpand θmare continuous variables,while the last onejis the discrete variable.If the discrete variable is reformulated into the innermost level of this multilevel optimal problem,the use of KKT condition will circumvent the complex derivation of constraints on independent variables (process parameters or model parameters).Compared with the nested reformulation and the decoupled reformulation,the computational amount of this method can be decreased intuitively.

If the number of process parameters ispand the number of model parameters isq,the total number of all potential uncertain parameters isp+q=k.In the traditional flexibility analysis,the uncertain set is ak-dimensional hyperrectangle.For example,ifp=1,q=1 thenk=2，as Fig.4(a) shows,the uncertain set is a defined rectangle.Whileviathe simplified reformulation,this 2-dimensional optimal problem can be simplified on the 1-dimensional number axis.Similarly,ifk=3，in Fig.4(b),the uncertain set is a 3-dimensional cuboid in θm-θp1-θp2coordinate system.If we use simplified reformulation,the θm-dimension will be eliminated and the optimal solution will only be found in the 2-dimensional plane.Fig.4(c) geometrically shows the dimension reduction of the simplified reformulation.

Fig.4.The simplified reformulations for developed flexibility analysis (a) example 1, p=1, q=1, k=2,(b) example 2, p=2, q=1, k=3 and (c) general situation.

2.3.Comparisons between the modified and existing reformulations

In order to assess the efficiency of the modified reformulation and to compare it with two existing reformulations,investigations on three cases are carried out.All cases are implemented in GAMS 24.9.2 and solved by global optimization solver ANTIGONE 1.1 on an Intel?CoreTMi5-8250U machine (2.30 GHz,8 threads,8 GB of RAM).Unless specified otherwise,problems are solved in parallel using all available threads.

2.3.1.Numerical case

The constraints of the numerical case are expressed as three linear inequalities,Eq.(20)to Eq.(22).The variation range of process parameter and model parameter are described by Eq.(23) and Eq.(24).

In Fig.5(a),inequality constraints,nominal point,and expected deviations are presented clearly.As shown in Fig.5(b),we expand the available operating range only in the θpdirection until one of constraints becomes active.Through the modified reformulation method conducted in Section 2.2,we can get the corresponding flexibility index δp=1.75 that is the same as the result of the work by Ochoaet al.[37].

Solutions from nested reformulation,recoupled reformulation,and the modified reformulation for this numerical case are shown in Table 1.

2.3.2.Deprotection reaction case

Deprotection reaction is a very common process in pharmaceutical industry.Fig.6 briefly shows the deprotection reaction in a continuous stirred tank reactor (CSTR).

The protecting group is cleaved from A to produce the desired product P and a gaseous byproduct B.It is assumed that a single phase within the reactor can be ensured by controlling the reaction pressure.Two main constraints of this process exist.Eq.(25) indicates that the conversion rate of A should be greater than 98%,while Eq.(26) emphasizes that the output concentration of P should be maintained greater than 0.45 mol·L-1.

Moreover,we can write down the mass balance equation and kinetic equation,where residence time τ is 350 min and temperature is 80 °C.

In this example,the process parameter refers to the input concentration of A,noted as.The permitted variation range ofis expressed by Eq.(32).Model parameters are pre-exponential factork0,and activation energyE,which are expressed by Eq.(33)and Eq.(34),respectively.The expectancy of both uncertain process and model parameters and their standard deviation are shown in Table 2.

Again,nested reformulation,recoupled reformulation,and the modified reformulation are applied.The corresponding results are shown in Table 3.The flexibility indexes come to consistence in δp=4.927.

Table 1Solutions of different reformulations for numerical case

Table 2The nominal values and the expected deviations for uncertain process and model parameters for the deprotection reaction case

2.3.3.Cascade reaction case

This case considers the cascade reactions in two CSTRs [38],as Fig.7 shows.

Process constraints include the minimum yield of product D and the minimum ratio of product D to unreacted ingredients.

Equalities regarding the kinetic and mass balance issues are represented by Eq.(37) to Eq.(39).In Eq.(37) and Eq.(38),r1andr2represent reaction rates.In Eq.(39),andcirepresent the input and output concentration of materiali.τ is the residence time.vijmeans the stoichiometric number of componentsiin reactionj.

Fig.5.The sketch map of numerical case.

Fig.6.The sketch map of deprotection reaction.

Uncertain model parameters are rate constantskj.The qualified expected range of τ andRB/Acan be expressed in Eq.(40) and Eq.(41),where Δτ is 90 min and ΔRB/Ais 1.δmensures model parameters remaining in a region with a cumulative probability that agrees with the desired confidence levelpc=0.8 [39].

Solutions are shown in Table 4.The flexibility index is δp=1.387,which is also as same as the work by Ochoaet al.[13].

Table 3Solutions of different reformulations for the deprotection reaction case

Table 4Solutions of different reformulations for the cascade reaction case

Table 5The nominal values and the expected deviations for uncertain process and model parameters for continuous ibuprofen manufacturing case

Fig.7.The sketch map of cascade reaction.

Benefits from the modified reformulation have been evaluated through three classical cases.Clearly,the modified reformulation can provide the same solution as those from the nested reformulation and recoupled reformulation but with less computational burden.Indeed,the size of the modified reformulation regarding the number of equations as well as the number of single and discrete variables is smaller.In the modified reformulation,the number of the discrete variables equals to the number of inequality constraints.The scale of equations and variables of the modified reformulations is about one third of that of nested reformulation.Such computational results are basically consistent with the theoretical derivations in Section 2.2.

3.Flexibility Analysis for Continuous Ibuprofen Manufacturing

The detailed continuous flow synthesis of ibuprofen is depicted in Fig.8[7].It contains three reactions including Friedel-Crafts acylation,1,2-aryl migration,and saponification reaction.Isobutylbenzene(IBB)and propionic acid are firstly mixed before the adding of triflic acid (TfOH).The mixtureF0is heated to 150 °C for meeting the requirements of the reaction happened in plug flow reactor(PFR)1.The outletF1,along with the mixture of rimethyl orthoformate (TMOF),MeOH and PhI(OAc)2serve as the feedstock for the reaction in PFR 2.Then MeOH and H2O are mixed to dissolve KOH and the solution is mixed withF3to formF4in PFR 3.The final step is to saponify methyl ester with KOH at the proper temperature 65 °C to get the crude productF5.Ibuprofen aqueous phase is finally generated after the purification ofF5.The detailed model,parameters and scalars of the continuous ibuprofen manufacturing are provided in Supplementary Material.

3.1.Model formulation for the flexibility analysis

Eq.(44) enforces the inequality constraints.g1andg5are in charge of ensuring the quality of production,whileg2,g3andg4guarantee the safety requirements.

Fig.8.The continuous-f low synthesis of ibuprofen.

By bringing the detailed equations from Supplementary Material Eq.(S30) to Eq.(S63) into Eq.(45),w e have the follow ing inequality equations.

Recalling Section 2.1,it is clear that the above f ive inequality constraints constitutegj（d，z，x，θ）in the flexibility analysis model,w hile all equations of thermodynamics and kinetics of reactions,mass balance,and heat balance given in appendix B Eq.(S30) to Eq.(S63) constitute thehi（d，z，x，θ）.In addition,state variablesxcan be eliminated as the principle that Eq.(3)and Eq.(4)illustrate.Consequently,gj（d，z，x，θ）andhi（d，z，x，θ）can be merged intof j（d，z，θ）.As an example,ifFA，0is chosen as the control variable,the corresponding flexibility analysis model can be formulated as Eq.(46).Of course,the flexibility analysis model w ith other choices of control variables can be formulated similarly.

3.2.Flexibility analysis with the modified reformulation

In order to illustrate the modified reformulation that was discussed in Section 2.2,this subsection concentrates on the flexibility analysis for ibuprofen manufacturing process without any control variables.In the continuous ibuprofen manufacturing process,the process parameters are the temperature of reactor 1T1,the input concentration of IBBCA，0,the temperature of reactor 2T2and the input concentration of PhI(OAc)2CD，2.Two model parameters are the pre-exponential factorA1,and activation energyE1.The nominal point/expectancy and expected/standard deviation of uncertain parameters are shown in Table 5.

Table 6The nominal values and variations of uncertain process parameters

Hence,these qualified intervals of process parameters can be expressed through Eq.(47) to Eq.(52).The considered model parameters are presented in Eq.(51) to Eq.(52) as it is assumed that δm=1.

The modified reformulation constrained by above inequalities for the flexibility analysis can be expressed by Eq.(53).The inequality constraints are the same as that in Eq.(46).Eq.(47) to Eq.(52) constitute the constraints of both process and model uncertain parameters.As no control variables are considered in this model,we pre-defineFA，0andFD，2are respectively 50 mmol·h-1and 38.865 mmol·h-1which are the same values as scenario 2 in the next subsection.

It can be easily observed that all reformulations again come to a consistence at δp=0.367.Furthermore,the total CPU time associated with the modified reformulation can be reduced by 1-3 orders of magnitude when compared to those with the nested reformulation or the recoupled reformulation.It should be noted that if the nested reformulation method is used,a number of bilinear terms will appear and will lead to a great computational cost.Indeed,for some cases,the optimization can’t be even converged within the pre-defined 1000 s.

3.3.Influences of different combinations on control variables

In order to simplify the complexity of the problem,following cases consider only uncertain process parameters while taking control variables into account.For the investigated continuous ibuprofen synthesis process,there are four potential uncertain process parameters -the temperature of reactor 1T1,the input concentration of IBBCA，0,the temperature of reactor 2T2,and the input concentration of PhI(OAc)2CD，2.Their variations are given in Table 6.Alternative control variables are the mole flow of IBBFA，0,and the mole flow of PhI(OAc)2FD，2.Apparently,they can be adjusted to cope with changes on uncertain parameters.Five scenarios are examined,respectively,-without any control variables,only with the control variableFA，0,only with the control variableFD，2,and with bothFA，0andFD，2.

Results are shown in detail in Table 7.Note that under the first two scenarios,no control variable means thatFA，0andFD，2need to be pre-defined.Otherwise,they are variables that need to be optimized.For scenarios with no control variable,the flexibility index is 0.411 and 0.786,respectively.In these two scenarios,no recourse action can be taken so that the flexible region is limited and the flexibility index is low.Interestingly,compared with the values of flexibility index listed in Table 8,scenario 1 and scenario 2 in Table 7 have larger flexibility index.The reason is that Section 3.2 considered both uncertain process parameters and model parameters,i.e.,constraints on uncertain model parametersA1andE1were requested.While Section 3.3 only considers uncertain process parameters.

Table 7The results of different combination scenarios of flexibility analysis models in ibuprofen case

Table 8Solution results and solution reports of different methods in continuous ibuprofen manufacturing case

In scenario 3,whenFA，0is considered as the control variable andFD，2is pre-defined with 38.865 mmol·h-1,FA，0can be optimized with 48.444 mmol·h-1.In this scenario,the flexibility index can be increased to 1.234 compared with scenario 2 because of the consideration of recourse actions.For scenario 4,ifFA，0is fixed as same as scenario 1,FD，2can be solved to optimality.Scenario 4 tells us,benefited from the recourse actions,the flexibility index can be increased to 0.867.In scenario 5,FA，0andFD，2are considered as the control variables simultaneously,in other words,recourse actions can be more flexible.As a result,the flexibility index can reach a high value at 2.104.

Taking scenario 4 as an example,Fig.9 reveals the relations between the flexibility index and control variables.GivenFA，0,FD，2can be correspondingly optimized to hedge against effects caused by uncertain process parameters.Different value of fixedFA，0and adjustedFD，2can provide corresponding distinct flexibility properties with different indexes.Fig.9 shows the trends of regulating effect of control variableFD，2,from which it indicates that theadjustment ofFD，2makes a great difference on flexibility under the variety ofFA，0.Needless to say,the relations illustrated by Fig.9 can help decision makers optimize the operation conditions under uncertainty so that the flexibility property can be guaranteed.For example,the decision makers can determine the value ofFA，0according to the expected flexibility index.

Fig.9.The flexibility index in different set FA，0 with control variable FD，2.(a) and (b) are from different angels of view.

However,it should be noted here that the results discussed here do not indicate that decision makers must select or select more control variables.Whether the control variables are considered or not depends on the real practical production needs.For decision makers,if the process system requires high flexibility,control variables should be selected without any doubt.However,if there is no strong requirement for flexibility which means the inherent flexibility has already met the production needs,the control variables won’t need to be considered.

4.Conclusions

This paper concentrated on the flexibility analysis of the continuous ibuprofen manufacturing process.Based on the existing nested reformulation and decoupled formulation,a modified reformulation for flexibility analysis was discussed to moderate the computational burden.Three classical cases were carried out to successfully illustrate the effectiveness of the modified reformulation.In detail,the model sizes and the CPU times were less than those of the existing nested reformulation and decoupled formulation.For the investigated continuous ibuprofen manufacturing process,CPU time associated with the modified reformulation was reduced by three orders of magnitude.Parallel to identifying the flexibility index without any consideration of control variables,different combinations on control variables were investigated to reveal the benefits from taking recourse actions into account.In general,such relations between the flexibility index and choice of control variables set a solid basis for further optimizing operation conditions.

In the future,in addition to the introduction of control variables,uncertain process and model parameters will be simultaneously considered to evaluate the flexibility property of continuous pharmaceutical manufacturing.Furthermore,the limitation of the modified reformulation is its precondition of independent uncertain parameters.Elliptical uncertain set instead of only hyperrectangle uncertain set will be investigated to reveal the effects of coupled uncertainties.Clearly,mathematical formulations will be updated and more efficient solution methods will be accordingly developed.On the other hand,much more uncertainties will also be considered in our future work.However,the more uncertain parameters are considered,the greater computational time will be cost.Needless to say,the tailored algorithm will be proposed.

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.

Acknowledgements

The authors gratefully acknowledge the financial support from the National Key Research and Development Program of China(2018AAA0101602).

Supplementary Material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cjche.2021.10.019.