Dynamic modelling and optimisation of the batch enzymatic synthesis of amoxicillin

: Amoxicillin belongs to the β-lactam family of antibiotics, a class of highly consumed pharmaceutical products used for the treatment of respiratory and urinary tract infections, and is listed as a World Health Organisation (WHO) “Essential Medicine”. The demonstrated batch enzymatic synthesis of amoxicillin is composed of a desired synthesis and two undesired hydrolysis reactions of the main substrate (6-aminopenicillanic acid (6-APA)) and amoxicillin. Dynamic simulation and optimisation can be used to establish optimal control policies to attain target product specification objectives for bioprocesses. This work performed dynamic modelling, simulation and optimisation of the batch enzymatic synthesis of amoxicillin. First, kinetic parameter regression at different operating temperatures was performed, followed by Arrhenius parameter estimation to allow for non-isothermal modelling of the reaction network. Dynamic simulations were implemented to understand the behaviour of the design space, followed by the formulation and solution of a dynamic non-isothermal optimisation problem subject to various product specification constraints. Optimal reactor temperature (control) and species concentration (state) trajectories are presented for batch enzymatic amoxicillin synthesis. dynamic optimisation


Introduction
Antibiotics are societally essential pharmaceutical products, making previously untreatable illnesses such as pneumonia and tuberculosis curable, thus revolutionising modern medicine [1]. Access to essential medicines in low-to middle-income countries and antibiotic shortages in developed countries remain an important issue [2,3]. The complex molecular structures of antibiotics imply expensive, multistep and materially intensive syntheses required to make such molecules [4]. Designing efficient and cost effective antibiotic manufacturing routes is imperative [5].
The family of β-lactam antibiotics includes some of the most important pharmaceutical products, with cephalosporins and semisynthetic penicillins corresponding to around 65% of the global production of antibiotics [6]. Figure 1 illustrates the leading antibiotics by share of infection treatments in the U.K. in 2016, with five of the leading nine antibiotics belonging to the β-lactam family. Table 1 lists certain β-lactam antibiotics with their applications, average sales volumes [7,8] and unit prices, taken from the National Chemical Database Service. The production of high-salevolume β-lactam antibiotics is typically implemented via enzymatic methods [9]. Amoxicillin is one β-lactam antibiotic listed as a World Health Organisation (WHO) "Essential Medicine": It is applicable to a variety of ailments [10], including respiratory and urinary tract infections (UTIs) and is the top antibiotic in terms of dosage worldwide [11]. The demonstrated enzymatic synthesis of amoxicillin paves the way for the elucidation of optimal design and operating parameters via modelling and optimisation [12,13].
Dynamic modelling and optimisation of batch processes have been substantially implemented in the literature for a variety of bioprocesses, including fermentations [14,15], antibiotic synthesis [16,17] and many other applications, illustrating the utility and versatility of such methods in bioprocess design. Enzymatic synthesis, crystallisation and reactive crystallisation studies for βlactam antibiotics have been implemented in the literature [16][17][18][19]. A demonstrated batch enzymatic synthesis of amoxicillin paves the way for theoretical studies [20]. Dynamic modelling and optimisation of the batch enzymatic synthesis of amoxicillin has yet to be implemented. Moreover, introducing temperature dependence into the model may allow for the optimisation of batch reactor temperature profiles to meet a specific production objective, which, to the best of our knowledge, has not been previously implemented. The main objectives and novelty of this work are as follows: 1. To introduce temperature dependency into the published kinetic model of batch enzymatic amoxicillin synthesis; 2. To understand the attainable performances and inherent trade-offs (via isothermal operation) of varying batch times and operating temperatures; 3. To optimise dynamic temperature profiles toward optimal process performance for varying product quality constraints. First, the enzymatic synthesis pathway is presented with the kinetic model. The temperature dependence of kinetic parameters is introduced through the regression of Arrhenius constants from published experimental data. An isothermal simulation is then implemented for design space investigation. A constrained dynamic optimisation problem is then described, which compares the resulting temperature (control) profiles required to meet certain product specification objectives pertinent to the specific process.

Amoxicillin Synthesis Pathway and Kinetic Model
The kinetically controlled reaction pathway for the synthesis of amoxicillin catalysed by penicillin G acylase (PGA, denoted as E) is shown in Figure 2. A classical mechanistic approach to the complete mechanism for this system has been shown to likely lead to an intractable problem: Thus, a simplified semi-empirical kinetic model is used. Amoxicillin (AMOX) is synthesised from the reaction of p-hydroxyphenyl glycine methyl ester (PHPGME) and 6-aminopenicillanic acid (6-APA). There are two side reactions: PHPGME hydrolysis to p-hydroxyphenyl glycine (PHPG) and amoxicillin hydrolysis to 6-APA and PHPG. Pertinent species properties are listed in Table 2. The enzyme concentration, CE, equals that considered in the literature [20].  Broad n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a Equations (1)-(8) describe the kinetic model for the reaction pathways shown in Figure 2 [20]. Here, C = the species concentration at time t; vi = the species rate of formation; vh1 = the rate of PHPGME hydrolysis; vh2 = the rate of AMOX hydrolysis; vS = the rate of AMOX synthesis; k = the species inhibition constant; kEN = the 6-APA adsorption constant; kcat = the reaction rate constant; KM = the reaction empirical rate constant; and Xmax is the maximum conversion ratio of the enzyme reagent complex into AMOX. Subscripts i and j denote species and reactions, respectively. The system of dynamic ODEs is solved simultaneously using the built-in MATLAB ODE solver ode15s.

Kinetic Parameter Estimation
In this section, temperature dependency is introduced to the published kinetic model through the regression of Arrhenius parameters for certain parameters in Equations (1)- (8). The considered operating conditions of the batch reactor for enzymatic amoxicillin synthesis are as described in the experimental demonstration in the literature [22]. A 2:1 molar mixture of PHPGME and 6-APA (40:20 mM, respectively) in a 1-L batch volume was used for the batch enzymatic synthesis, with no amoxicillin or PHPG present at the start of the reaction, i.e., CPHPG (t0 = 0), CAMOX (t0 = 0) = 0. The original dynamic model for the batch enzymatic synthesis of amoxicillin assumes isothermal operation, i.e., it does not account for the temperature dependence of kinetic parameters. Here, the kinetic parameter temperature dependence for amoxicillin is introduced from published temperature-dependent data: Values of kcat,1 and kcat,2 were regressed from published amoxicillin concentration data at 5, 25, and 35 °C [22], from which Arrhenius parameters could then be estimated. The remaining kinetic parameters were assumed to be temperature-independent, as a much wider kinetic dataset is required for further multiparametric regression. Values of kcat,1 and kcat,2 at different temperatures were regressed by minimising the residual error between the kinetic model (Equations (1)-(8)) and the experimental data using the bound constrained solver "fminsearchbnd" in MATLAB. Table 3 shows the regressed values of kcat,1 and kcat,2 at the given temperatures. With temperature-varying values for kcat,1 and kcat,2, the Arrhenius parameters, k0 (the preexponential factor) and Ea (the energy barrier), were then regressed. An Arrhenius-type temperature dependence of kcat,1 and kcat,2 was assumed, according to Equation (9): where k0 and Ea are the Arrhenius pre-exponential factor and energy barrier, respectively; R is the universal gas constant; and T is the reaction temperature. The fitting methodology for Arrhenius constant regression also used "fminsearchbnd" in MATLAB, as described previously. Figure 3 shows the lines of best fit for both Arrhenius plots, showing a good fit in both cases. The regressed Arrhenius parameter values (listed in Table 4) allowed for good replication of the experimental data ( Figure 4). A corroboration of the regressed parameters with a wider dataset will further validate the values used in this work. The experimental data in the literature that were used to regress the kinetic parameters at different operating temperatures (T = 5, 25, 35 °C) provided no error values for the calculated concentrations (estimated via HPLC) or temperatures, and thus error bars are not shown [22]. An investigation of the effect of errors on regressed parameter values (isothermal kcat,1 and kcat,2 and non-isothermal Arrhenius parameters) could be implemented by perturbing their values and observing the effect on the optimisation results as a form of sensitivity analysis.

Design Space Investigation and Simulation
Surface plots of the concentrations of key species, PHPGME, 6-APA, AMOX and PHPG as a function of isothermal reactor temperature and time were developed for design space investigation. Additionally, various performance indices were investigated and compared to elucidate attainable performances and inherent trade-offs. Selectivity, S, is the ratio of concentration of AMOX to PHPG (Equation (10)), and productivity, P, is the rate of production of AMOX (Equation (11)). Surface plots of these performance indices were also developed to understand the behaviour of the dynamic system.

Dynamic Optimisation
The optimisation problem is posed such that different target amoxicillin product concentrations are met while consumption of the valuable feedstock 6-APA (bulk price = 35-40 USD kg −1 ) is minimised. Minimising feedstock consumption in different cases can lead to possible process configurations where reagents and solvent are recycled following desired product (amoxicillin) crystallisation [23,24]. Although recycling options are not explicitly considered in this work, they have been demonstrated in the literature for β-lactam antibiotic synthesis with nanofiltration for byproduct (PHPG) removal [25]. This can allow for processes with greater material efficiencies and overall sustainability. Economic analyses of specific process designs can further elucidate the benefits of such designs [5], but are beyond the scope of this work. The problem considers system state variables, i.e., species concentration profiles (C(t)) that vary over time and are influenced by the reactor temperature control vector T(t). The problem aims to simultaneously consider the minimisation of 6-APA consumption with the scope of potential recycling of this key substrate along with amoxicillin product maximisation. The bicriteria problem is defined by Equations (13) and (14). While J1 and J2 present competing objectives, the utility of the presented optimisation results is in exploring how to manipulate the optimal temperature control profile over the batch duration to yield the desired amount of amoxicillin while minimising feedstock consumption. Equations (12)-(14) are Numerous approaches can be used to modify a multiobjective problem for compatibility with single-objective solution methods. Commonly, a weighted sum objective is used to combine competing objectives into a single term, with weights defining the relative importance of each. However, the weights assigned to various process targets to produce a single objective function may be considered arbitrary in many cases, with decision-makers not necessarily able to quantify a priori the relative importance of competing objectives. Rather, we elected to consider an ε-constraint approach [26]. One of the objectives can be considered to be a constraint in the problem formulation, Model (35 C) while the other is solved to optimality. This is repeated, increasing the value of the objective constraint by ∆ε and re-solving, repeating this process by incrementally increasing the constraint value across the entire span of permissible values for that particular objective. In this case, the amoxicillin product concentration was treated as a secondary objective and converted into a constraint. The single objective herein was to minimise the consumption of key feedstock 6-APA (formulated as the maximisation of final 6-APA concentration), subject to varying final amoxicillin concentrations (end-point constraints) and a final batch time, tf = 500 min, which was observed as the maximum time beyond which amoxicillin degradation via hydrolysis to 6-APA and PHPG dominated (Figures 3 and 4). Equations (12)- (14) become Equations (15)- (18): s.t.
Reactor temperature control trajectories were considered to be piecewise constant with N = 20 time discretisation elements. A variety of different initialisation profiles were tested to investigate the sensitivity of the problem to the initial guess of temperature trajectory, which consistently converged to the same optima presented in this work. Although global optima were not guaranteed, a variety of initialisation profiles resulting in the same optima and profiles implied that as close to global optima were attained as possible with the solution methods implemented. Temperature control profiles were initialised as isothermal at T = 295 K across all time discretisations. Lower and upper boundaries on temperature values were 5 and 35 °C, respectively (interior constraints, Equation (18)), corresponding to the temperature range used to regress kinetic parameter Arrhenius constants. Different end-point constraints on amoxicillin concentration (Equation (16)) were considered, as well as an unconstrained case (ε = 0). Each case of a considered amoxicillin concentration constraint was solved as a separate optimisation problem. The problem was solved using a direct, simultaneous method: Orthogonal collocation on finite elements [27,28] to approximate the state (C(t)) and control (T(t)) trajectories to convert the problem into a nonlinear programme (NLP), which was then solved using an interior point filter line search algorithm (IPOPT) [29]. Dynamic optimisation was implemented in MATLAB using DynOpt [30].

Dynamic Simulation
Concentration surface plots of PHPGME, 6-APA, AMOX and PHPG as a function of batch time and reactor operating temperature are shown in Figure 5. The colours in Figure 5 represent concentration data values on the z axis and are shown to aid in the interpretation of the surface plots. Species PHPGME and PHPG concentrations were strong functions of batch runtime but were largely unaffected by temperature, which was likely due to the model being fitted exclusively to AMOX concentration data [22]. Key substrate (6-APA and PHPGME) concentrations were lowest at longer batch times and higher temperatures, as these conditions favoured consumption toward amoxicillin and its hydrolysis to PHPG. Amoxicillin synthesis was favoured by longer batch times and higher reactor temperatures, while PHPG concentration was not affected significantly by reactor temperature. The concentration of AMOX reached a maximum at ~430 min, after which it began to decrease due to the decreasing ratios of vS to vh2 and vAMOX to vPHPG (as shown in Figure 6), both of which decreased significantly with batch time. The colours in Figure 6 represent the reaction rate ratios shown on the z axis and are shown to aid in the interpretation of the surface plots.

Non-Isothermal Simulation
Extensive simulations were performed as a preliminary investigation of attainable process performances for non-isothermal batch reactor operations. Here, the T(t) domain was discretised on a coarse grid (discretisation level N = 9) to generate a finite set of profiles for exhaustive simulation. Possible temperature profiles were subject to the constraint that the temperature was allowed to change by a maximum of 10 °C per interval. The different simulation cases are presented in Figure 7, showing the maximum attainable amoxicillin concentrations versus the maximum batch time required (tMAX), and also selectivity (Equation (10)) and productivity (Equation (11)). The banding observed was due to stepwise temperature profile simulations. Different cases are highlighted on each of these plots: (1) maximum selectivity, (2) maximum amoxicillin concentration and (3) a compromise between selectivity and productivity. Maximum final amoxicillin concentration was achieved by operating at the maximum temperature (35 °C), in agreement with previously observed results [22]. There are inherent trade-offs between different process performance metrics. Doing so facilitates a visualisation of the attainable performance of the process, subject to the rules imposed in generating the set of profiles. Dynamic optimisation can be implemented for temperature profile manipulation in order to investigate the possibility of the process benefitting from non-isothermal operation for specific production objectives.

Non-Isothermal Dynamic Optimisation
The resulting temperature (control) and species concentration (state) trajectories for different product amoxicillin constraints imposed on the dynamic optimisation problem are shown in Figures  8 and 9, respectively. The temperature profiles varied significantly as the constraint concentration of amoxicillin in the reactor product mixture increased. When the optimisation problem was unconstrained, the temperature gradually decreased from 296 K at t = 0 min to 278.15 K at t = 325 min (lower temperature bound), after which isothermal operation at this temperature continued for the remainder of the batch time. For CAMOX (tf) = 5 and 6 mM, the temperature profile decreased until t = 375 min and then increased again, with a higher final operating temperature required for a higher concentration constraint. When the constraint was set to CAMOX (tf) = 7 mM, a gradual decrease from 296 to 290 K over the batch duration was implemented. For CAMOX (tf) ≥ 8 mM, temperature profiles increased over the batch duration. At CAMOX (tf) = 9 mM, a steady increase from 295 to 305 K was observed. For CAMOX (tf) = 10 mM, the temperature increased sharply from 292 to 308 K (upper temperature bound) at t = 100 min. For CAMOX (tf) = 11 mM, the temperature increase rate was more drastic, increasing to the upper bound at t = 25 min, with subsequent isothermal operation for the remainder of the batch duration. These trends confirmed that lower amoxicillin target constraints required lower operating temperatures, with higher temperatures required for higher target concentrations. This behaviour was due to the effect of increasing temperatures favouring amoxicillin synthesis over undesired side hydrolysis, consistent with the system behaviour shown earlier in Figure 5. For the unconstrained case, the temperature was eventually pushed to the lower bound, and as the target amoxicillin concentration increased, the temperature was pushed to the upper bound. The temperature bounds imposed on the dynamic optimisation problem corresponding to the lower and upper temperature values at which experimental concentration data for kinetic and Arrhenius parameter regression were implemented [22]. Relaxing these bounds may have resulted in different optimal control trajectories. Industrial temperature controllers typically operate similarly to piecewise constant or linear profile behaviours, but there is inevitably a temperature lag between the constant temperature discretisations. The incorporation of time lag into temperature control will certainly affect the presented results: Such an analysis was outside the scope of this work, but may enhance our understanding of attainable optima for this work and other industrially relevant design cases. The optimal temperature profiles for each case presented are those attained from our optimisation, having initialised the solver with the same starting solution (isothermal at T = 295 K). There were different profiles that could still meet the same amoxicillin product concentration constraints (i.e., different routes to the same product specification). The isothermal concentration response surfaces in Figure 5 were unable to depict state evolution under non-isothermal reactor operation, and thus dynamic optimisation was implemented.
The resulting species concentration (state) trajectories are shown in Figure 9. In all cases, as substrates PHPGME and 6-APA were consumed, product amoxicillin and by-product PHPG were formed. For the unconstrained and lower amoxicillin concentration constraints (e.g., CAMOX (tf) = 5-7 mM), the target amoxicillin concentration was met before the concentration of 6-APA (i.e., the objective function) was maximised. This indicated that there was scope to also optimise the batch runtime in the dynamic optimisation problem as well, which can be considered in future work. The final species concentrations at the end of the batch runtime (tf = 500 min) are shown in Figure 10. As the product amoxicillin constraint was increased, the maximum objective function value decreased. The objective of the dynamic optimisation was to minimise 6-APA consumption (posed by maximising the final 6-APA concentration at the end of the batch duration). Figure 10 shows the fraction of 6-APA from the process feed that was present in the batch product versus that consumed in the reaction network. In all design cases (varying the amoxicillin product concentration constraint), a significant portion of the fed 6-APA was present in the product mixture. It can be seen from Figure 9 that when CAMOX (tf) = 11 mM was specified, the constraint was not met. Consideration of the maximum attainable objective function values for a broad range of CAMOX (tf) (Figure 11) explains this result. This was a result of there being a maximum product concentration of amoxicillin attainable from the given feed substrate concentrations. Figure 12 plots the attained maximised objective function for different imposed amoxicillin product concentration constraints. Until a certain amoxicillin product concentration (Point A), the maximum attainable objective    The initial concentration values of PHPGME and 6-APA used in this work were the same as those implemented in the literature, from which experimental data were used for parameter regression. These values corresponded to typical industrial application ranges [20]. In the literature, it was observed that concentrations of 6-APA > 20 mM led to a slower formation of the acyl-enzyme complex in the amoxicillin synthesis reaction, which was detrimental to process performance [20,22]. For this reason, the initial conditions (i.e., reagent concentrations) considered in this work were kept as demonstrated in the literature. Additionally, an investigation of the applicability of the model at higher reactor volumes for increased production scales would be of value, requiring further experimental validation and corroboration. The results presented in this work assumed perfect mixing and heat transfer in the batch reactor. As the process is scaled up, deviations from this ideal behaviour assumption must be considered when interpreting the results presented in this work.
Measuring the material efficiency of different design cases is also an important consideration in pharmaceutical production [31]. A variety of green chemistry metrics exist for quantification of the effectiveness of processes, the applicability of which depends on the process [32]. Here, we compare the efficiency of different optimisation scenarios (varying CAMOX (tf) considerations) via the reaction mass efficiency (RME), calculated by Equation (19).
The RME is calculated as the mass ratio of amoxicillin (mAMOX) at the end of a batch run to the masses of starting materials (PHPGME (mPHPGME) and 6-APA (m6-APA)) at the start (t0 = 0). Values of RME for different product amoxicillin concentration constraints are shown in Figure 13. As the specified amoxicillin concentration constraint increases, the RME increases, which is expected. All values of RME were relatively midrange with respect to typical pharmaceutical manufacturing processes [33]. The effect of scale-up on material efficiencies, as well as plant-wide costs, is an important consideration during process design and optimisation studies such as this one [34].  Figure 14 shows purity, with different species as the reference point corresponding to optima (maximum product 6-APA concentrations) for different product amoxicillin concentration constraints, calculated from the final concentration values shown in Figure 10. In all cases, PHPGME and 6-APA (reagent) concentrations decreased as amoxicillin purities increased. While PHPG contents decreased, they remained high compared to other species in the mixture. A consideration of mixture compositions following synthesis is important for purification and separation process design, particularly for crystallisation design [35]. Understanding the partitioning of different species between liquid phase mixtures and solid phase product crystals is essential for understanding the effects of different operating and design parameters on crystal product attributes (purity, size distribution, etc.). Another method of circumventing the accumulation of impurities in crystals is reactive crystallisation, which aims to preferentially crystallise the desired compound from the reaction mixture solution: This has been implemented for various β-lactam antibiotics [16,17]. While such methods may allow for improved purities over traditional crystallisation, ensuring rates of 14 reaction compared to mass transfer rates to avoid high supersaturation is paramount, as this can lead to undesirably wide crystal size distributions and low mean sizes. Figure 14. Product mixture content for different amoxicillin product concentration constraints (colour scheme congruent with Figures 9 and 11).
The basis of the dynamic optimisation cases considered in this work was to establish optimal temperature profiles to meet specific amoxicillin product concentrations while minimising the consumption of key feedstock (6-APA). Varying constraints on the amoxicillin product concentrations were considered to investigate the effect of different desired mixture qualities on optimal temperature control trajectories. The objective of the dynamic optimisation problem was not to maximise amoxicillin production, but rather to produce at least a certain desired quantity of amoxicillin while minimising consumption of the key reagent 6-APA. In the presented results, dynamic temperature profiles were confirmed to be optimal for attaining a range of target amoxicillin concentrations while minimising 6-APA consumption. The highest target amoxicillin concentration considered had a nearly isothermal temperature profile over the batch duration: If the objective of the optimisation had been simply the maximisation of amoxicillin product concentration, isothermal profiles would likely have been preferable. However, for certain desired amoxicillin product concentrations, dynamic temperature profiles were confirmed to be optimal.

Conclusions
A kinetic model for the batch enzymatic synthesis of amoxicillin was considered in this work for dynamic modelling and optimisation. The novelty of this work is in the following: 1. The introduction of temperature dependency to the implemented existing kinetic model [20] using experimental concentration data [22] for the batch enzymatic synthesis of amoxicillin; and 2. The first implementation of dynamic temperature profile optimisation in order to meet specific product quality constraints and minimise feedstock consumption in β-lactam antibiotic production. Temperature dependency in the model was introduced through the regression of Arrhenius parameters to allow for non-isothermal modelling and optimisation of the temperature trajectories to minimise key substrate (6-APA) consumption subject to varying constraints on the product amoxicillin concentration. The attainable minimal substrate consumption was strongly dependent on the amoxicillin product concentration. This work presents the first batch reactor temperature