Next Article in Journal
Toward Optimized 89Zr-Immuno-PET: Side-by-Side Comparison of [89Zr]Zr-DFO-, [89Zr]Zr-3,4,3-(LI-1,2-HOPO)- and [89Zr]Zr-DFO*-Cetuximab for Tumor Imaging: Which Chelator Is the Most Suitable?
Next Article in Special Issue
Solid Lipid Microparticles by Spray Congealing of Water/Oil Emulsion: An Effective/Versatile Loading Strategy for a Highly Soluble Drug
Previous Article in Journal
Locust Bean Gum Nano-Based Hydrogel for Vaginal Delivery of Diphenyl Diselenide in the Treatment of Trichomoniasis: Formulation Characterization and In Vitro Biological Evaluation
Previous Article in Special Issue
Complexation: An Interesting Pathway for Combining Two APIs at the Solid State
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Pharma 4.0-Artificially Intelligent Digital Twins for Solidified Nanosuspensions

by
Christina Davidopoulou
1,2 and
Andreas Ouranidis
1,2,*
1
Department of Chemical Engineering, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece
2
Department of Pharmaceutical Technology, School of Pharmacy, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece
*
Author to whom correspondence should be addressed.
Pharmaceutics 2022, 14(10), 2113; https://doi.org/10.3390/pharmaceutics14102113
Submission received: 8 September 2022 / Revised: 28 September 2022 / Accepted: 30 September 2022 / Published: 3 October 2022
(This article belongs to the Collection Feature Papers in Pharmaceutical Technology)

Abstract

:
Digital twins capacitate the industry 4.0 paradigm by predicting and optimizing the performance of physical assets of interest, mirroring a realistic in-silico representation of their functional behaviour. Although advanced digital twins set forth disrupting opportunities by delineating the in-service product and the related process dynamic performance, they have yet to be adopted by the pharma sector. The latter, currently struggles more than ever before to improve solubility of BCS II i.e., hard-to-dissolve active pharmaceutical ingredients by micronization and subsequent stabilization. Herein we construct and functionally validate the first artificially intelligent digital twin thread, capable of describing the course of manufacturing of such solidified nanosuspensions given a defined lifecycle starting point and predict and optimize the relevant process outcomes. To this end, we referenced experimental data as the sampling source, which we then augmented via pattern recognition utilizing neural network propagations. The zeta-dynamic potential metrics of the nanosuspensions were correlated to the interfacial Gibbs energy, while the density and heat capacity of the material system was calculated via the Saft-γ-Mie statistical fluid theory. The curated data was then fused to physical and empirical laws to choose the appropriate theory and numeric description, respectively, before being polished by tuning the critical parameters to achieve the best fit with reality.

Graphical Abstract

1. Introduction

A digital twin is the mirrored functional space model of a manufacturing unit that responds to physical state changes which is purposed for surveillance, optimization and prediction [1].The depth of the replication appointment to the physical asset by the digital model space remains ambiguous, hence various integration levels occur, whilst optimality appears controversially dependent of the corresponding functionality, complexity and availability of infrastructures [2]. The birth of the digital twin pertained to a crafty inspiration of the NASA Apollo 13 mission controllers 52 years ago. In a need-based, lifesaving attempt, the crew set up and modified simulations to mimic the multiplicity of the spacecraft’s physical conditions, those which occurred 45,000 miles from the Earth’s surface. This first successful approach spurred inspiration across applied scientific fields.
Since then, digital modeling is emancipating the state-of-the-art approach for delivering process and product lifecycle awareness, enabling unparalleled plantwide control, optimization, and prediction for material manufacturing [3]. Akin to such product material ontologies, artificial neural networks (ANNs) have ascended as a surrogate, responsive framework, poised to identify and simulate non-linear dependencies between digital twin variables [4]. As such, ANNs constitute the requirements for demanding, multiple, complex experiments towards the generation of the product manufacturing cycle (PML) data, irrelevant. Moreover, ANNs recognize the relationship patterns between independent and dependent variables of the PML, whether the available data set is insufficient or noisy [5]. The latter advancement proves to be of extreme importance, hence such variables lie at the core of the digital twin mechanistic and/or empirical algorithms, whilst poor data set availability effects common burdens when dealing with expensive raw and starting materials, such as those handled by the pharma industry. Successful applications of the ANNs in the field include the calculation of API’s critical quality attributes and their physicochemical properties, the prediction of the in-vitro drug release profile, the identification of raw material–tablet properties characterization, the dissolution behavior, and the size effect of injectable microparticle prediction [6].
For such a heavily regulated sector, ANN-fused artificially intelligent digital twins might accurately predict the critical process parameters, the final product’s quality attributes, integrate the process steps and lead seamless scale-up studies. The advantages of possessing this blend of stochastic, empirical, and mechanistic knowledge becomes invaluable in terms of technical, economical and risk eliminating factors by the creation of novel, in-silico intelligent analogues to the physical product. Counterintuitively, although digital twin applications could transform the pharmaceutical industry, whereas precision and risk elimination are evidently required, the penetration of these Pharma 4.0 tools appears limited. In addition, artificially intelligent, advanced digital twin systems, i.e., ANNs-fused digital twins, have not been reported by the field’s literature.
Under this lens, arguably, one significant contemporary technological challenge of the pharmaceutical industry is to improve solubility and in vivo dissolution profiles of poorly soluble active pharmaceutical ingredients (APIs). To address this challenge, APIs are mixed with stabilizers and formulation excipients composites and comminuted by wet media milling into nanosuspensions. The obtained, coated, liquid crystals are further stabilized by spray drying solidification to allow for further processability. Taking into account the proven reliability and scalability of this method, we introduced novel thermodynamic models to calculate such composite material solubility by assessing the obtained stabilizer-coated nanoparticle Gibbs energy anisotropic minimization, quantified by the implementation of PC-SAFT interrogations coupled with elastic tensor analysis [7].
In this current research, we build on these foundations to create a novel, fully operational, neural network-intelligent digital twin, capable of describing the course of manufacturing of solidified nanosuspensions given a defined PML launching chronic point, and also predict and optimize the engaged process outcomes [2]. Elaborating on this strategy, the ANN was embedded in a stepwise manner among data sampling, model deployment and curve fitting towards the implementation of the digital twin (Figure 1). The ANN fulfilled the mission of augmenting the abundance of the available discrete data generated to calibrate and validate the model whilst eliminating the experimental burden and the model’s uncertainty [8].

2. Materials and Methods

2.1. Process Design

The API’s formulation processes, including wet milling and spray drying steps, were modelled using the Siemens gProms Formulated Products® platform Process Systems Enterprise, gPROMS, (www.psenterprise.com/products/gproms, accessed on 2 August 2022). Before formulation, it is crucial to select the compatible polymer or surfactant to coat the liquid crystals, the one which best enhances the API’s dissolution performance [9]. This pivotal step was included in the process model; the selection was based on the implementation of the computational statistical associating fluid theory criteria (PC-SAFT) to substitute for the experimental trials.
In detail, the stabilizer’s addition generates a semi-solid interface that adds to the dissolution’s Gibbs energy decrease by GEE (J mol−1), according to Equations (1)–(3) [7].
G E E = G m s + G m i = R T l n K 2 K 1
G m s = 2 γ V m r ( 1 C r )
G m i = 1.7 ε A P I σ A P I ρ s t a b Δ ( σ s t a b A P I ) ( ε s t a b A P I ) m s t a b γ V m r
The G m s and G m i terms refer to the particle’s surface and the interface caused by the stabilizer Gibbs energy. Where K 2 K 1 is the ratio of the dissolution equilibrium coefficient post and prior the powder’s size decrease respectively, T (K) is the dissolution temperature and R (J mol−1 K−1) is the universal gas constant. In Equation (2) γ (Ν m−1) is the surface tension, Vm (m3 mol−1) the API’s molar volume, r (m) the particle’s characteristic size and C (m) a parameter equal to 1.5 ( V m N A ) 1 3 , while in Equation (3) Δ (m) is the material’s distance between the molecular layers, mstab (−) the stabilizer’s number of segments per chain based on PC-SAFT theory, and ε (eV) and σ (m) are the depth of pair potential and the segment diameter, respectively. The terms εstab-API and σstab-API were calculated using the Berthelot-Lorentz combining rules. Stabilizer candidates were chosen, namely Poloxamer-188, Poloxamer-407 and HPC-SL. As far as poloxamers are concerned, since both are copolymers containing ethylene oxide (EO) and propylene oxide (PO) groups in different proportions, the component’s PC-SAFT parameters were calculated using Berthelot-Lorentz rules (see Table 1).
The polymer that contributes dominantly to the Gibbs energy decrease is the one to be chosen as the stabilizer. Gibbs energy is not calculated experimentally, yet on the other hand the zeta-dynamic potential can be correlated to the interfacial Gibbs energy. The physical parameter estimations regarding the API and the stabilizers were performed using the gProperties® package of Siemens Process Systems Enterprise (https://www.psenterprise.com/products/gproms/properties, accessed on 2 August 2022). The Saft-γ-Mie equations of state were utilized to calculate the temperature dependence on the density and the heat capacity of the components.
The wet milling process model was based on the main batch grinding Equation (4) and three supplementary empirical grinding functions, each of them complemented by their corresponding breakage rate. This batch grinding mass balance set up explains the change of the powder’s mass fraction containing particles of size interval i:
d w i d t = j = 1 i 1 [ S j w j d ( b ( i , j ) ) d x i ] S i w i  
The initially examined breakage function was the one proposed by de Vegt et al. [15], and it is referred to as the de Vegt model (Equation (5)). The corresponding breakage rate is strongly dependent on the product’s material properties (Equation (6)):
b ( i , j ) = S i S j
S ( i ) = c   E k i n , i E f r a c t , i P y ρ V H x i K 1 c
E k i n , i = W m , k i n ρ V i
E f r a c t , i = 0.896 ( π ( 1 υ 2 ) Υ ) 2 3 ( 0.0183 δ 2 ( V 0 V i   ) 1 / 4 ) 5 / 3
In Equations (5)–(7), S(i) (s−1) is the breakage rate of a particle of size interval i, Ekin (J) and Efract (J m−3) are the kinetic energy of the particles and the fracture energy, respectively, Py (Pa) is the yield pressure, ρ (kg m−3) is the particle’s density, V (m3) is the mill’s chamber volume, H (Pa) is the particle’s hardness, xi (m) is the particle size i, K1C (Pa m−1/2) is the stress intensity factor, Wm,kin the mass specific impact energy (J kg−1) and b(i,j) (−) is the mass fraction of the product that fell from size interval j to i. In Equation (8), δ is the solubility parameter (Pa1/2), Vi (m3) is the particle’s volume, υ (−) is the Poisson’s ratio, Υ (Pa) is the Young’s modulus of elasticity and V0 (m3) is the unit’s crystal volume. While all the other parameters in Equations (6) to (8) are pre-estimated, the breakage rate parameter c (−) is a tuning parameter, estimated by experimental data. The second breakage function and breakage rate used for the wet milling simulations is proposed by Austin et al. [16] and it is referred to as the Austin model (see Equations (9) and (10)):
b ( i , j ) = φ ( x i x j ) γ + ( 1 φ ) ( x i x j ) β
S ( i ) = {         a ( x i x c r i t ) d       x j x c r i t                     0                 x j < x c r i t  
where φ, γ, β, d (−) and a (s−1) are the tuning parameters, while xi and xj are the product’s and the post-breakage final particle size accordingly (m), and xcrit is the critical particle size, namely the size after which no breakage occurs. The last breakage function was the one proposed by Kapur et al. [17] and it is referred to as the Kapur model (model Equations (11) and (12)):
b i , j = ( x i x j ) e
S i = A x i k
where e, k (i) and A (s−1) are the tuning parameters. Apparently, all three breakage functions and rates encumber tuning and physical parameters respectively, with Austin presenting the highest number of considered tuning parameters and de Vegt the lowest. The parameter estimation was conducted in the Siemens Process Systems Enterprise gPromsFP® Model Validation platform (https://www.psenterprise.com/products/gproms/modelbuilder, accessed on 2 August 2022), applying the Maximum Likelihood Estimation method. In the spray drying process, the droplets’ hydrodynamic diameter was assumed to obey lognormal distribution and for the drying rate calculation the Oakley’s model was adapted [18]. For each particle size interval i was defined from the particle size distribution and was dispersed in a droplet, which in turn belongs to a size interval z, the local mass and energy balance describing the spray drying model which is described respectively by the Equations (13) and (14).
m ˙ s , i , z d x i , z , j , t d t = N ˙ i , z , j , t
m ˙ s , i , z ( C p , s , i + j = 1 N j x i , z , j , t C p , j ) d T i , z d t = h A i , z ( T g T i , z ) + m ˙ s , i , z j = 1 N j λ j d x i , z , j , t d t  
where m ˙ s , i , z (kg s−1) is the corresponding solids particles flowrate, x i , z , j , t (kg kg−1) is the dry basis moisture content of the liquid specie j, N ˙ i , z , j , t (kg s−2) is the drying rate time derivative, C p , s , i (J kg−1 K−1) is the solid material’s specific heat capacity and C p , j the liquid specie’s corresponding one, T i , z (K) is the solid particle’s temperature, h (J m−2 s−1 K−1) is the heat transfer coefficient, A i , z (m2 s−1) is the shrinking rate of the surface area of the droplet, T i , z (K) is the droplet’s temperature, and λj (J kg−1) is the latent heat of vaporization of the liquid specie j. The local vapor phase’s mass balance for the evaporated liquid specie j and the vapor’s phase energy balance is described in Equations (15) and (16), respectively.
d m v , j d t = m ˙ v , i n x v , j , i n m ˙ v , o u t x v , j , o u t + i = 1 N i z = 1 N z t = 0 t τ N ˙ i , z , j , t d t
d H v d t = H v , i n H v , o u t + i = 1 N i z = 1 N z j = 1 N j t = 0 t τ N ˙ i , z , j , t d t   C p , v , j ( T v T 0 ) i = 1 N i z = 1 N z h t = 0 t τ A i , z , t ( T v T i , z ) d t
where m ˙ v , i n and m ˙ v , o u t (kg s−1) are the inlet and outlet mass flowrate of the vapor phase respectively, x v , j , i n and x v , j , o u t (kg kg−1) are the inlet and outlet mass fractions of the liquid specie j in the vapor phase respectively, H v , i n and H v , o u t (J s−1) are the inlet and outlet enthalpy flowrates of the vapor phase respectively and C p , v , j (J kg−1 K−1) is the specific heat capacity of liquid specie j in the vapor phase, tτ (s) is the droplets residence time inside the spray dryer’s chamber and Tv (K) is the vapor phase’s temperature. The unhindered drying rate Νu,z,j (kg s−1), i.e., the drying rate of the very same droplets without containing solids described by Equation (17) and is interlinked to the actual drying rate by the relative drying rate f (−) (Equation (18)).
N ˙ u , z , j , t = h λ j A i , z , t ( T v T w b , j )
N ˙ i , z , j , t = f N ˙ u , z , j , t
where Twb,j (K) is the wet bulb temperature of the liquid specie j.

2.2. Experimental Study and Digital Twin Thread Structuring

The wet milling process of the API Itraconazole nanosuspension was performed by a Pulverisette 7 Premium (Fritsch GmbH, Idar-Oberstein, Germany). Delivered discrete experimental values [19] were first used to train an ANN and form a complete size reduction profile (see Table 2), while afterwards this profile is used as a basis to fit the tuning parameters of the digital twin model. The process reduced the powder’s particles’ sizes approximately by one order [8]. Post wet milling, the micronized API suspension was spray dried to remove the liquid phase and obtain the desired dry powder product [19]. The spray dryer used in the experiment was the Büchi B-191 Mini Spray dryer (Büchi, Flawil, Switzerland).
The technical specifications of the mill’s mechanical parameters, such as the rotation and the revolution speed, the capacity and the equipment’s volume, were determined by the manufacturer. Moreover, the specific ranges of input data were adopted by the manufacturer’s technical specification sheet of the Buchi B191 Mini Spray dryer (Büchi, Switzerland). Both data input sets are herein demonstrated in Table 3.
The system descriptive equations were concluded by Section (B), containing linear and differential equations. For instance, the Equations (13)–(18) were considered as a system of six non-linear algebraic equations:
f ( t ,   x ˙ i , z , j , t ,   N ˙ i , z , j , t ) = 0      Equation (13)
f ( t , x i , z , j , t , x ˙ i , z , j , t , T i , z , T ˙ i , z ) = 0      Equation (14)
f ( t , m ˙ v , j , N i , z , j , t ) = 0      Equation (15)
f ( t , H ˙ v , N i , z , j , t , T i , z ) = 0      Equation (16)
f ( t , N ˙ u , z , j ) = 0      Equation (17)
f ( t , N ˙ u , z , j ,   N ˙ i , z , j , t ) = 0      Equation (18)
The system’s variables embedded first order derivatives and were calculated utilizing finite difference approximations. Specifically, the Backward Differentiation Approximations (BDF) were adopted, which present a universally approved method to approach DAE system solutions [20]. By applying the BDF using a time step size hs, at time tn, the DAE system transforms into a linear system featuring six equations bearing six unknown variables, each derivative of them being approximated as in Equation (25).
  x ˙ i , z , j , t n = h s 1 l = 0 k a l x i , z , j , t n l
where k is the degree of the interpolating polynomial and al is the lth polynomial’s level coefficient.

2.3. Integration of Artificial Neural Networks for Parameter Tuning

The ANN was properly trained to simulate the dynamic comminution profile inside the physical mill. Data sampling during the wet mill process is a difficult task, hence it requires temporary terminations. In addition, the final crystal size itself is not adequate for the characterization of the mill’s performance. It is therefore crucial to obtain additional information of the dynamic profile in order to fulfil the purposed optimization and prediction purposes. The discrete experimental training data includes the D50 values taken after sampling in between 6-min intervals. The transfer function used for the hidden layers was the Sigmoid Function (Equation (26)), due to its ability to identify non-linear relationships. For the output layer, the transfer function selected was the Rectified Linear Unit (ReLU), hence the case was addressed through a regression scenario and not as a classification (Equation (27)). The ANN training adopted the error backpropagation methodology towards defining the weights’ adjustments, but its layout was forward propagating. Furthermore, the weights wi,j connecting a layer with nj neurons with the next layer with ni neurons, followed the He initialization (Equation (28)).
f ( x ) = 1 1 e x
g ( x ) = m a x   ( 0 , x )
w i , j   2 n j           i ( 1 , n i )   ,     j ( 1 , n j )
For the identification of the D50(t) profile, polynomial regression was used with high efficiency, especially whereas polynomial fitting is of a higher order (≥4).

3. Results

3.1. Material Critical Quality Attributes and the Material System’s Interfacial Gibbs Energy Assesment

The density and the heat capacity temperature profile of the pure stabilizers, calculated via the Saft-γ-Mie statistical fluid theory, are shown in Figure 2. Poloxamer-188 material poses the lowest density and the lowest heat capacity, while HPC is denser than poloxamers. Calculations of the pure stabilizers’ densities were pivotal for the determination of the stabilizer’s interfacial Gibbs energy. Apparently, density effects the intermolecular forces of the stabilizer-API composite (see Table 4). For higher temperatures, the composite’s cohesive forces decrease as the kinetic energy of the molecules increases, and thus they tend to escape their structured positions [21]. For the same reason, when temperature increases the material tends to expand, the density decreases, and the corresponding particle’s surface tension γ (Ν m−1) decreases as well, see Equation (29) [9]. In addition, the zeta-potential bourn within the semi-solid interface that the stabilizer’s addition forms [22], is utilized as indicator of its efficiency (Table 3). A stabilizer causing high absolute values of zeta-potential, creates strong repulsive forces (Coulomb forces) preventing the particles from aggregating.
γ = 0.33 k B T ( N A ρ M r ) 2 3 [ ln ( S 0 55.6 ) + 5 ]
In Equation (29), where kB (J K−1) is the Boltzmann constant, Mr (g mol−1) and ρ (g m−3) are respectively the molecular weight and the density of the particle, NA (mol−1) the Avogadro number and S0 (−) the solubility of the pure API in water in absence of the stabilizer.
According to the results shown in Table 4, the free Gibbs energy of the interface is correlated both to the density and the zeta-potential. The decrease of surface density and tension enhances the binding with the stabilizer, which in turn plays a crucial role in the zeta-potential absolute value increase [23]. While Gibbs energy is a measurement of the maximum non-expansion work available in a system, interfaces of higher Gibbs free energy favors the effects of the electrostatic Coulomb forces. Among the three investigated stabilizers, Poloxamer-188 presents the highest interfacial Gibbs energy, and as a result the highest contribution in the dissolution Gibbs energy enhancement (Equation (1)) and the highest absolute zeta-potential, making it suitable for selection.

3.2. Parameter Fitting

The results of the milling’s dynamic profile mapping appear in Figure 3. The artificial neural network achieved fitting within 1.6% mean squared error (MSE), while the second-order polynomial achieved fitting with 1.1% MSE. Although using polynomial fitting provides a bit lower percentage of MSE, the advantages of using ANNs for pattern recognition are numerous, as discussed in the introduction section. Also, the discrete experimental data points appear to be noisy. Deep neural networks identify noisy data rather than memorize it and include it in the training process [24]. However, there exist data fusion algorithms with insignificant computational requirements used to minimize noise during sampling, such as the Extended Kalman filter.
The Maximum Likelihood Estimation best-fit results comparing the three breakage functions are shown in Figure 4. The Austin model provided the best experimental fit with the lowest MSE, while the de Vegt model provided the highest one. Furthermore, as discussed above, the Austin model’s breakage function included more tuning parameters than the rest, with the Kapur model being the second and the de Vegt model being the last one. Considering this fact, it was found that the polyparametricaly tuned model exhibited a suitable experimental fit curve. As expected, the integration of numerous tuning parameters enhances the curve fitting performance, especially when it desired to interpolate complex experimental data profiles. The tuning parameter herein successfully plays the role of “correction factor”.

3.3. Sensitivity Analysis

3.3.1. Wet Milling

A sensitivity analysis was performed against the various considered tuning parameters regarding the breakage functions conducted to examine their effects on the model outputs. Figure 5 shows the effect of the de Vegt model’s breakage rate parameter c (−) on the final D50 size profile. Parameter c is the tuning parameter existing in the model, while the others are predetermined according to the material’s and the mill’s characteristics (see Equations (6)–(8)).
The de Vegt model’s breakage distribution function remained the same during the analysis as when combining Equations (5) and (6) the final form appears dependent on the x i x j size reduction ratio (Equation (30)).
b ( i , j ) = ( x i x j ) 1.25
Apparently, increasing c provides higher breakage rates for larger particles and results in efficient comminutions. However, with the latter featuring the only tuning parameter for a given API, each desired final size shall unveil its own unique reduction profile. De Vegt model’s breakage function is ideal for qualitative analysis cases, e.g., when only the initial and the final sizes are of the main interest. In Figure 6 and Figure 7, Austin’s comminution profiles are illustrated, based on the breakage rate’s sensitivity analysis. The effects of the tuning parameters a (s−1) and d (−) were interrogated, and it was found that by increasing each or both of a and d, the breakage rates of interval sizes i (Equation (10)), and thus the comminution efficiency, increases (Figure 6). The breakage distribution function b(i,j) during the Austin model’s relative sensitivities remained the same, with φ = 0.3, γ = 1.17 and β = 4. In comparison to de Vegt model, a desired D50 size is achieved via various profile paths (Figure 7a). In addition, even when φ = 1 and γ = 1.25, namely when the distribution function b(i,j) is the same with the de Vegt’s one, the existence of the exponent d as a secondary tuning parameter in the S(i) calculation function allows the formation of multiple profiles leading to the same result (Figure 7b). This is the benefit of engaging multiple tuning parameters within the breakage rate function, as they render the digital twin capable of fitting multiple profiles.
Figure 8 and Figure 9 illustrate the Kapur model D50 time profile results (e = 6), following the same strategies. The Kapur model proved to be a sufficient breakage function as it recognizes multiple time profiles. Nevertheless, the b(i,j) distribution function contains one tuning parameter in comparison with the Austin model, and this is the reason that effective experimental fit accuracy cannot be approached by the theory.

3.3.2. Spray Drying

Two-factor sensitivity analysis was conducted to the spray drying model to map the design space. Figure 10 presents the obtained analysis’s contour diagram, whereas the factors considered are the air temperature and the air volume flow rate, while the response is the final product’s humidity, which constitutes the critical quality attribute related to the process. The contours below describe the final moisture content of the product decrease as the air’s temperature and/or its volumetric flow rate increases for the given initial dry humidity (g g−1). The air water capacity threshold rises proportionally with temperature when the temperature increases the air’s given humidity (g g−1) covering a lower proportion of the moisture threshold and its relative humidity decreases, allowing humidity absorption [7].

4. Discussion

In our artificially intelligent digital twin thread, comminution efficiency appears dependent both on the material system selection and on the equipment’s parameter settings. The relevant capacity and the physical properties of the API and of the grinding media were included in the digital twin’s working parameters, defining the quality attributes of the final product. Therefore, this digital twin thread is capacitated to analyze what-if equipment and material scenarios for early production assessments towards the reduction of material waste and the optimization of time schedule and process functions serving multiple batch mill units [25]. In addition, the algorithm is poised to identify irreversible faults and incompatibilities beforehand, for the related processes and the material system (e.g., over-efficient comminution profiles, API’s unwanted loss by dissolution during milling etc.) and propose controller actions to avoid them. Spray drying is a complex procedure, crucial for the final API’s powder formation. It is a continuous process enabling real-time data sampling and immediate response to changes of the physical object’s properties. In our digital twin, the spray drying process is described as a multi-parametric block, since the final temperature and moisture content of the product depend on the drying air’s temperature, pressure, and capacity as far as the initial product’s flow, temperature, and droplet size distribution. Therefore, it can predict the drying efficiency given the air temperature specification, hence apart from the convergent mass and energy balances it exploits psychrometric calculations to determine the air’s relative humidity and its moisture capacity threshold [7], thus being advantageous in delivering realistic analysis reports.
Our algorithm may receive real-time data and return forecasts of the temperature profile inside the spray dryer for any given initial input conditions, an extremely important operation, hence the product is required to be obtained in a stable crystalline state, while its temperature should not overcome the API’s and the excipients’ melting point [26,27]. This digital twin ameliorates any scale-up related efforts and product degradation risks, while contributing to the optimization of the process outlining the system’s design spaces. In this recursive digital model, the proper empirical law selection, the one that more accurately describes the physical object’s progress, is dependent on the crowd of its tuning parameters, namely the parameters acceptable to be adjusted within boundless values towards improvement of the model’s prediction accuracy. Considering processes like wet milling, although the material and energy balances invoke fundamental laws of physics, models that describe specific functions, such as the comminution laws, must include data-driven parameters. Under this lens, for this specific unit block, using trained neural networks could replace such models efficiently, although limitations should be considered, and those lie in the idiosyncratic nature of the methods such as the ability of the empirical models to point out relative dependencies between variables [28]. ANNs identify patterns without obeying a strict mathematical function, thus generalizing the relationship between independent-dependent variables. Both ANNs and polynomial fitting are efficient, depending on the circumstances chosen. In addition, as the adjective “empirical” implies, experiments cannot be avoided completely, hence even the structure of the mathematical models requires real data for their parameter training.
Contemplating the limitations of the integration depth of our approach, there exist cases where a specific level of uncertainty is required, for instance when model discrepancy needs to become recognized. A model can simultaneously be mechanistically biased, including over-confident parameter estimates, and therefore can also be effected by model discrepancy. In this research, the adjusted parameters were not predetermined, hence they do not possess physical status, yet their existence lies within calibration purposes, such as fitting experimental data towards the enhancement of the digital twin’s validity. Such tuning parameters, although as iterated irrelevant of physical interrelations, become scientifically vital, whereas complex process simulations are considered, as they aid the translation of the data set ontologies to the related physical property [8]. Adjusting tuning parameters to discrete data points vertically amplifies the model’s data value and therefore improves overall accuracy as uncertainty, and its minimization, shall be in frame with the output results.
The digital twin thread requires input modifications when APIs other than Itraconazole are considered. Specifically, different APIs provide different experimental data and consequently different ANN weights and physichochemical parameters must be plugged to the system’s agent. ANNs, however, carry the flexibility to link different input types with desired outputs, and as such various physicochemical properties can be utilized, constituting this digital twin as a universal predictive platform suitable for BCS II drugs [29], whereas the iterated processes are examined. Pharma 4.0 capacitating digital twins fused with artificial wisdom will in the future decrease the need for the implementation of experiments and exploit real time manufacturing process data towards the constant adjustment of the tuning parameters, succeeding the control and prediction of the process’s progress and material’s quality attributes in time. Although fully functional digital twins require demanding computational and sensoring infrastructure as continuous curve fitting, probabilistic forecast, information evaluation and system configuration are taking place simultaneously, their multimodal capabilities shall be undeniably useful for industry deployment in the years to come.

Author Contributions

Conceptualization, A.O. and C.D.; methodology, A.O. and C.D.; software, A.O. and C.D.; validation, A.O. and C.D.; formal analysis, A.O. and C.D.; investigation, A.O. and C.D.; resources, A.O. and C.D.; data curation, A.O. and C.D.; writing—original draft preparation, A.O. and C.D.; writing—review and editing, A.O. and C.D.; visualization, A.O. and C.D.; supervision, A.O.; project administration, A.O.; funding acquisition, A.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The authors wish to acknowledge the Siemens Process Systems Engineering company and support team for the knowledge-based advice.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

  1. Kritzinger, W.; Karner, M.; Traar, G.; Henjes, J.; Sihn, W. Digital Twin in Manufacturing: A Categorical Literature Review and Classification. IFAC-PapersOnLine 2018, 51, 1016–1022. [Google Scholar] [CrossRef]
  2. Udugama, I.A.; Lopez, P.C.; Gargalo, C.L.; Li, X.; Bayer, C.; Gernaey, K.V. Digital Twin in Biomanufacturing: Challenges and Opportunities towards Its Implementation. Syst. Microbiol. Biomanuf. 2021, 1, 257–274. [Google Scholar] [CrossRef]
  3. Ouranidis, A.; Davidopoulou, C.; Tashi, R.K.; Kachrimanis, K. Pharma 4.0 Continuous MRNA Drug Products Manufacturing. Pharmaceutics 2021, 13, 1371. [Google Scholar] [CrossRef] [PubMed]
  4. Sutariya, V.B.; Groshev, A.; Pathak, Y.V. Artificial Neural Networks in Pharmaceutical Research, Drug Delivery and Pharmacy Curriculum. In Proceedings of the 29th Southern Biomedical Engineering Conference 2013, Miami, FL, USA, 3–5 May 2013; pp. 91–92. [Google Scholar] [CrossRef]
  5. Mandlik, V.; Bejugam, P.R.; Singh, S. Application of Artificial Neural Networks in Modern Drug Discovery; Elsevier Inc.: Amsterdam, The Netherlands, 2016; ISBN 9780128015599. [Google Scholar]
  6. Wang, S.; Di, J.; Wang, D.; Dai, X.; Hua, Y.; Gao, X.; Zheng, A.; Gao, J. State-of-the-Art Review of Artificial Neural Networks to Predict, Characterize and Optimize Pharmaceutical Formulation. Pharmaceutics 2022, 14, 183. [Google Scholar] [CrossRef]
  7. Ouranidis, A.; Davidopoulou, C.; Kachrimanis, K. Integrating Elastic Tensor and Pc-Saft Modeling with Systems-Based Pharma 4.0 Simulation, to Predict Process Operations and Product Specifications of Ternary Nanocrystalline Suspensions. Pharmaceutics 2021, 13, 1771. [Google Scholar] [CrossRef]
  8. Brynjarsdóttir, J.; Ohagan, A. Learning about Physical Parameters: The Importance of Model Discrepancy. Inverse Probl. 2014, 30, 114007. [Google Scholar] [CrossRef]
  9. Ge, K.; Ji, Y.; Lu, X. A Novel Interfacial Thermodynamic Model for Predicting Solubility of Nanoparticles Coated by Stabilizers. Chin. J. Chem. Eng. 2021, 31, 103–112. [Google Scholar] [CrossRef]
  10. Ferreira-Pinto, L.; De Araujo, P.C.C.; Aranda Saldaña, M.D.; Arce, P.F.; Cardozo-Filho, L. Experimental Data and Thermodynamics Modeling (PC-SAFT EoS) of the {CO 2 + Acetone + Pluronic F-127} System at High Pressures. J. Chem. Eng. Data 2019, 64, 2186–2192. [Google Scholar] [CrossRef]
  11. Bodratti, A.M.; Alexandridis, P. Formulation of Poloxamers for Drug Delivery. J. Funct. Biomater. 2018, 9, 11. [Google Scholar] [CrossRef]
  12. Li, B.; Zhang, L.; Zhang, Z.; Gao, R.; Li, H.; Dong, Z.; Wang, Q.; Zhou, Q.; Wang, Y. Physiologically Stable F127-GO Supramolecular Hydrogel with Sustained Drug Release Characteristic for Chemotherapy and Photothermal Therapy. RSC Adv. 2018, 8, 1693–1699. [Google Scholar] [CrossRef] [Green Version]
  13. Szafraniec, J.; Antosik, A.; Knapik-Kowalczuk, J.; Chmiel, K.; Kurek, M.; Gawlak, K.; Odrobińska, J.; Paluch, M.; Jachowicz, R. The Self-Assembly Phenomenon of Poloxamers and Its Effect on the Dissolution of a Poorly Soluble Drug from Solid Dispersions Obtained by Solvent Methods. Pharmaceutics 2019, 11, 130. [Google Scholar] [CrossRef] [Green Version]
  14. Guth, F.; Staal, B.; Jackson, T. Kolliphor ® P 188 Bio: A New Era in Shear Protection; BASF Corporation: Florham Park, NJ, USA, 2017; pp. 1–4. [Google Scholar]
  15. De Vegt, O.; Vromans, H.; Faassen, F.; Van Der Voort Maarschalk, K. Milling of Organic Solids in a Jet Mill. Part 2: Checking the Validity of the Predicted Rate of Breakage Function. Part. Part. Syst. Charact. 2006, 22, 261–267. [Google Scholar] [CrossRef]
  16. Austin, L.; Shoji, K.; Bhatia, V.; Jindal, V.; Savage, K.; Klimpel, R. Some Results on the Description of Size Reduction as a Rate Process in Various Mills. Ind. Eng. Chem. Process Des. Dev. 1976, 15, 187–196. [Google Scholar] [CrossRef]
  17. Kapur, P.C. Self-Preserving Size Spectra of Comminuted Particles. Chem. Eng. Sci. 1972, 27, 425–431. [Google Scholar] [CrossRef]
  18. Oakley, D.E. Spray Dryer Modeling in Theory and Practice. Dry. Technol. 2004, 22, 1371–1402. [Google Scholar] [CrossRef]
  19. Pavlidou, M.; Κyriakos, K. Influence of Formulation Factors on Production of Microcomposite Particles for Inhalation; Aristotle University of Thessaloniki: Thessaloniki, Greece, 2018. [Google Scholar] [CrossRef]
  20. Brenan, B.K.E.; Engquist, B.E. Backward Differentiation Approximations of Nonlinear Differential/Algebraic Systems; American Mathematical Society: Providence, RI, USA, 1998; Volume 51, pp. 659–676. Available online: Https://www.Jstor.Org/Stable/2008768 (accessed on 8 September 2022).
  21. Thanuja, B.; Kanagam, C.; Sreedevi, S. Studies on Intermolecular Interaction on Binary Mixtures of Methyl Orange-Water System: Excess Molar Functions of Ultrasonic Parameters at Different Concentrations and at Different Temperatures. Ultrason. Sonochem. 2011, 18, 1274–1278. [Google Scholar] [CrossRef]
  22. Kojima, T.; Karashima, M.; Yamamoto, K.; Ikeda, Y. Combination of NMR Methods to Reveal the Interfacial Structure of a Pharmaceutical Nanocrystal and Nanococrystal in the Suspended State. Mol. Pharm. 2018, 15, 3901–3908. [Google Scholar] [CrossRef]
  23. Shnoudeh, A.J.; Hamad, I.; Abdo, R.W.; Qadumii, L.; Jaber, A.Y.; Surchi, H.S.; Alkelany, S.Z. Synthesis, Characterization, and Applications of Metal Nanoparticles; Elsevier Inc.: Amsterdam, The Netherlands, 2019; ISBN 9780128144282. [Google Scholar]
  24. Rolnick, D.; Veit, A.; Belongie, S.; Shavit, N. Deep Learning Is Robust to Massive Label Noise. arXiv 2017, arXiv:1705.10694. [Google Scholar]
  25. News, T.; News, T. Improving Pulp M Ill Ill Operation Operation With Digital Twin Twin Technology. Tech News 2018, 2, 49–50. [Google Scholar]
  26. Ouranidis, A.; Gkampelis, N.; Markopoulou, C.; Nikolakakis, I.; Kachrimanis, K. Development of a Nanocrystal Formulation of a Low Melting Point Api Following a Quality by Design Approach. Processes 2021, 9, 954. [Google Scholar] [CrossRef]
  27. Singh, A.; Van den Mooter, G. Spray Drying Formulation of Amorphous Solid Dispersions. Adv. Drug Deliv. Rev. 2016, 100, 27–50. [Google Scholar] [CrossRef] [PubMed]
  28. Scardi, M. Artificial Neural Networks as Empirical Models for Estimating Phytoplankton Production. Mar. Ecol. Prog. Ser. 1996, 139, 289–299. [Google Scholar] [CrossRef] [Green Version]
  29. Ouranidis, A.; Tsiaxerli, A.; Vardaka, E.; Markopoulou, C.K.; Zacharis, C.K.; Nicolaou, I.; Hatzichristou, D.; Haidich, A.-B.; Kostomitsopoulos, N.; Kachrimanis, K. Sildenafil 4.0—Integrated Synthetic Chemistry, Formulation and Analytical Strategies Effecting Immense Therapeutic and Societal Impact in the Fourth Industrial Era. Pharmaceuticals 2021, 14, 365. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Digital twin deployment and parameter estimation strategy by steps: (a) real data sampling from the physical object, (b) data multiplication via pattern recognition using artificial intelligence algorithms aiming accuracy enhancement, (c) determination of the proper descriptive physical and empirical laws, and (d) adjustment of tuning parameters to achieve best curve fit.
Figure 1. Digital twin deployment and parameter estimation strategy by steps: (a) real data sampling from the physical object, (b) data multiplication via pattern recognition using artificial intelligence algorithms aiming accuracy enhancement, (c) determination of the proper descriptive physical and empirical laws, and (d) adjustment of tuning parameters to achieve best curve fit.
Pharmaceutics 14 02113 g001
Figure 2. Saft-γ-Mie stabilizer candidates, material properties results (a) Poloxamers’ densities (b) Poloxamers’ heat capacities (c) HPC-SL density (d) HPC-SL heat capacity as a function of temperature.
Figure 2. Saft-γ-Mie stabilizer candidates, material properties results (a) Poloxamers’ densities (b) Poloxamers’ heat capacities (c) HPC-SL density (d) HPC-SL heat capacity as a function of temperature.
Pharmaceutics 14 02113 g002
Figure 3. Wet milling dynamic comminution profile of experimental data points, ANN generated and best-fit polynomial curve.
Figure 3. Wet milling dynamic comminution profile of experimental data points, ANN generated and best-fit polynomial curve.
Pharmaceutics 14 02113 g003
Figure 4. Result dynamic profiles for the three breakage functions (a) de Vegt function’s best fit curve generating a 0.89% MSE (b) Kapur function’s best fit curve generating a 0.08% MSE, and (c) the Austin function’s best fit curve generating a 0.02% MSE.
Figure 4. Result dynamic profiles for the three breakage functions (a) de Vegt function’s best fit curve generating a 0.89% MSE (b) Kapur function’s best fit curve generating a 0.08% MSE, and (c) the Austin function’s best fit curve generating a 0.02% MSE.
Pharmaceutics 14 02113 g004
Figure 5. Breakage parameter c effect on the size comminution profile.
Figure 5. Breakage parameter c effect on the size comminution profile.
Pharmaceutics 14 02113 g005
Figure 6. D50 time profiles with Austin’s breakage function in various a and d parameter values with γ = 1.17 and φ = 0.3 (a) a = 8 × 10−5 and (b) d = 1.6.
Figure 6. D50 time profiles with Austin’s breakage function in various a and d parameter values with γ = 1.17 and φ = 0.3 (a) a = 8 × 10−5 and (b) d = 1.6.
Pharmaceutics 14 02113 g006
Figure 7. Austin’s breakage function for a D50 = 1 μm final size via various time profile paths with (a) γ = 1.17 and φ = 0.3 and (b) with γ = 1.25 and φ = 1 (de Vegt approximation).
Figure 7. Austin’s breakage function for a D50 = 1 μm final size via various time profile paths with (a) γ = 1.17 and φ = 0.3 and (b) with γ = 1.25 and φ = 1 (de Vegt approximation).
Pharmaceutics 14 02113 g007
Figure 8. D50 time profiles with Kapur’s breakage function in various k and A parameter values with e = 6 (a) A = 0.008 and (b) k = 1.2.
Figure 8. D50 time profiles with Kapur’s breakage function in various k and A parameter values with e = 6 (a) A = 0.008 and (b) k = 1.2.
Pharmaceutics 14 02113 g008
Figure 9. Kapur’s breakage function for a D50 = 1μm final size via (a) various time profile paths with e = 6 and (b) various profile paths with e = 1.25 (de Vegt approximation).
Figure 9. Kapur’s breakage function for a D50 = 1μm final size via (a) various time profile paths with e = 6 and (b) various profile paths with e = 1.25 (de Vegt approximation).
Pharmaceutics 14 02113 g009
Figure 10. Contour sensitivity analysis diagram of the spray drying model.
Figure 10. Contour sensitivity analysis diagram of the spray drying model.
Pharmaceutics 14 02113 g010
Table 1. PC-SAFT parameters of EO and PO groups for the calculation of the Poloxamers copolymers’ corresponding ones.
Table 1. PC-SAFT parameters of EO and PO groups for the calculation of the Poloxamers copolymers’ corresponding ones.
Groupmseg (−)σi (A)ui/k (K)Source
EO0.052 ΜWtotal2.89206.74[10,11,12,13,14]
PO0.037 MWtotal3.34192.72
Table 2. Discrete experimental sampling values used as training data for the ANN. These values represent three performed experiments using Itraconazole as API and Poloxamer-188 as stabilizer [19], which proved to be the most suitable after the Gibbs energy analysis.
Table 2. Discrete experimental sampling values used as training data for the ANN. These values represent three performed experiments using Itraconazole as API and Poloxamer-188 as stabilizer [19], which proved to be the most suitable after the Gibbs energy analysis.
Time (s)Experiment 1
D50 (μm)
Experiment 2
D50 (μm)
Experiment 3
D50 (μm)
3601.391.311.37
7200.8740.7920.846
14400.7840.7420.783
21600.5220.5010.520
28800.4670.4340.436
36000.3050.2970.301
Table 3. Input and output parameters considered for the development of the digital twin.
Table 3. Input and output parameters considered for the development of the digital twin.
MODELPARAMETERTYPEVALUEUNIT
Water quantityinput9mL
API contentinput0.5g
Stabilizer contentinput0.25g
Mannitol contentinput1g
Wet millInitial particle size (D50)input1.5μm
Grinding timeinput1h
Rotor speedinput600rpm
Rotor diameterinput40mm
Equipment volumeinput48mL
D50(t)output--
Air temperatureinput110°C
Air flowinput800L h−1
Spray dryerAir pressureinput5bar
Drying chamber volumeinput5L
Drying timeinput1h
Final product size (D50)output10μm
Table 4. Interfacial Gibbs energy and material properties correlation.
Table 4. Interfacial Gibbs energy and material properties correlation.
Stabilizer G m i Density (kg m−3)Z-Potential
HPC-SL0.00191320−11.7
Poloxamer-4070.0039954−13.7
Poloxame-1880.0056951−17.0
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Davidopoulou, C.; Ouranidis, A. Pharma 4.0-Artificially Intelligent Digital Twins for Solidified Nanosuspensions. Pharmaceutics 2022, 14, 2113. https://doi.org/10.3390/pharmaceutics14102113

AMA Style

Davidopoulou C, Ouranidis A. Pharma 4.0-Artificially Intelligent Digital Twins for Solidified Nanosuspensions. Pharmaceutics. 2022; 14(10):2113. https://doi.org/10.3390/pharmaceutics14102113

Chicago/Turabian Style

Davidopoulou, Christina, and Andreas Ouranidis. 2022. "Pharma 4.0-Artificially Intelligent Digital Twins for Solidified Nanosuspensions" Pharmaceutics 14, no. 10: 2113. https://doi.org/10.3390/pharmaceutics14102113

APA Style

Davidopoulou, C., & Ouranidis, A. (2022). Pharma 4.0-Artificially Intelligent Digital Twins for Solidified Nanosuspensions. Pharmaceutics, 14(10), 2113. https://doi.org/10.3390/pharmaceutics14102113

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop