A New Validation Methodology for In Silico Tools Based on X-ray Computed Tomography Images of Tablets and a Performance Analysis of One Tool

In silico tools which predict the dissolution of pharmaceutical dosage forms using virtual matrices can be validated with virtual matrices based on X-ray micro-computed tomography images of real pharmaceutical formulations. Final processed images of 3 different tablet batches were used to check the performance of the in silico tool F-CAD. The goal of this work was to prove the performance of the software by comparing the predicted dissolution profiles to the experimental ones and to check the feasibility and application of the validation concept for in silico tools. Both virtual matrices based on X-ray micro-computed tomography images and designed by the software itself were used. The resulting dissolution curves were compared regarding their similarity to the experimental curve. The kinetics were analysed with the Higuchi and Korsmeyers–Peppas plot. The whole validation concept as such was feasible and worked well. It was possible to identify prediction errors of the software F-CAD and issues with the virtual tablets designed within the software.


Introduction
In 2004, the American Food and Drug Administration (FDA) started the Process Analytical Technology initiative to transform the art of pharmaceutical development and processes to science with the goal of quality by design [1]. Since then the concept was revised several times. The International Council of Harmonisation adopted the concept for several guidelines [2,3]. Since the recent years, an increasing number of in silico tools to support the pharmaceutical development is available. However, a complete in silico development of pharmaceutical dosage forms does not seem possible at the present time. Nevertheless, in silico tools may contribute to cost reduction and saving of development time in the future.
Nowadays, different software packages like F-CAD (CINCAP, Basel, Switzerland) or DDDPlus (Simulations Plus Inc., Lancaster, CA, USA) provide assistance with the formulation development. The different in silico software tools are either based on cellular automata or the discrete element method and finite element method, respectively [4,5]. The example software F-CAD of this work used cellular automata and differential equations [6]. The concept of cellular automata was described by Wolfram (2002) and can be used to calculate complex correlations or phenomena with simple rule sets [7].
Virtual matrices were used as pharmaceutical dosage forms for the in silico calculations. The software can either design virtual dosage forms with a defined composition or upload Tagged Image File Format (TIFF) images as virtual matrix [6,[8][9][10]. Kimura et al. calculated the disintegration time of mefenamic acid containing tablets with a sufficient accuracy using F-CAD in 2013 [8]. The second goal of F-CAD is the prediction of dissolution profiles [6].
A reliable validation method for in silico tools is a mandatory step to enable a complete in silico development without experimental fitting for all pharmaceutical applications. The correctness of the prediction has to be verified. The opportunity of uploading TIFF images as virtual matrix for the in silico simulation provides the possibility to compare the predicted dissolution profile of X-ray micro-computed tomography (XµCT) images of dosage forms with the experimental dissolution profile of the imaged dosage forms. If both profiles match, the in silico tool can predict the behaviour of the API and the excipient(s) correctly. Furthermore, the virtual matrices designed by the software could be compared to the virtual matrices based on XµCT images to verify the correct virtual design of the tool.
If the prediction of a formulation is correct, the data can possibly be used to build up databases of excipients or placebo mixtures to transfer the prediction to new drugs, assuming negligible interactions between the different compounds. The potential was shown by Bollmann and Kleinebudde in a proof of concept study [11]. Due to limitations of the XµCT imaging contrasts between pharmaceutical powders, complex formulations could be accessible under the condition of negligible interactions between all excipients and using multiple binary mixtures for validation purposes. This would allow a complete in silico formulation development in line with the concept of quality by design [10,11]. However, the potential and the feasibility of this methodology needs to be further investigated.
The in silico tools require specific properties of XµCT images to upload them as virtual matrix. The image processing of the XµCT images has to guarantee that the resulting virtual matrix matches the real tablet. The pathway from the raw image to the uploadable TIFF image of tablets was presented for the tool F-CAD by Bollmann and Kleinebudde [9,10]. Some of the resulting images were used as virtual matrices for the prediction of dissolution profiles in this work [10]. The predicted dissolution profiles of the virtual matrices either based on the XµCT images or designed by the particle arrangement and compaction module (PAC-module) of F-CAD were compared to the experimental dissolution profiles.
Costa and Lobo reviewed several approaches to compare dissolution profiles [12]. In this work, the similarity ( f 1 and f 2 factor) of the predicted profiles compared to the experimental one is calculated as it is the recommended approach by the FDA and European Medicines Agency [12][13][14]. Furthermore, the kinetics are analysed with the Higuchi and Korsmeyer-Peppas plot to determine the release mechanism [15][16][17][18].
The goal of this work is to verify and rate the prediction of dissolution profiles of the software package F-CAD and to prove the concept and feasibility of the validation of in silico tools using XµCT images as virtual matrices. In addition, a sensitivity analysis with respect to virtual matrices based on XµCT images with different desirabilities is performed to prove the hypothesis of Bollmann and Kleinebudde about the impact of the desirability on the in silico dissolution calculation.

Excipients and Tableting
Theophylline anhydrate (BASF, Ludwigshafen, Germany, Dx50 114 µm) served as model active pharmaceutical ingredients (API) while ethyl cellulose (Aqualon N10, Ashland, Wilmington, DE, USA, Dx50 210 µm) was used as matrix former. Binary mixtures (100 g) in different ratios (25/75, 50/50, and 75/25) were blended with a Turbula mixer TC2 (Willy A. Bachofen AG, Muttenz, Switzerland) with 49 rpm for 15 min. Round, flatfaced tablets were compressed with a compaction pressure of approximately 160 MPa on a Styl'One Evolution (MEDEL'Pharm, Beynost, France). The target weight was 400 mg with a diameter of 11.28 mm. The tablets were already used in previous publications [9,10]. The data were used for additional analysis methods in this paper. The batches were named TxxEyy and the tablets were numbered. T25E75-01 means the tablet 1 of the batch containing 25 % theophylline and 75 % ethyl cellulose.

Characterisation
The helium density of both compounds was determined at 25°C using an AccuPyc 1330 (Micromeritics, Unterschleißheim, Germany (n = 3)). The tablets were characterised in height, diameter, and weight using a SmartTest50 (Sotax AG, Basel, Switzerland). More details are available in section 2.5 of the previous publication [9].
The solubility of the API was determined UV-spectroscopically (UV-1800 SHIMADZU, Kyoto, Japan). A suspension of theophylline in phosphate buffer pH 6.6 of the United States Pharmacopeia (USP) was stored and stirred in a water bath (37°C) for three weeks (n = 3). The suspension was filtered (0.2 µm, polypropylene), diluted and measured at a wavelength of 271 nm.

Dissolution Testing
Dissolution tests were performed with an automatic dissolution tester (DT 720, ER-WEKA, Langen, Germany) using apparatus 2 (USP) and 1000 mL phosphate buffer pH 6.6 (USP). The stirring speed was set at 100 rpm and the temperature was kept constant at 37°C. Samples were taken every ten minutes through frits (Poroplast, 10 µm). The drug release was determined UV-spectroscopically (UV-1800 SHIMADZU, Kyoto, Japan) in flow-through cuvettes at a wavelength of 271 nm. The total API content was determined by the plateau after complete dissolution of the tablets [9].

Imaging and Image Processing
The tablets were imaged by a CT-ALPHA (ProCon X-ray, Sarstedt, Germany). The amperage was 100 µA and the voltage 86 kV. In total, 1600 projections were taken while the sample rotated 360°with a voxel size of 7 µm.

Nomenclature and Methods of Pathways
The following notations were used for the different image processing pathways: bleach correction: reference (Ref) [20], stack (Sta) [21,22]; lighting correction: remove background (RB) [23], pseudo flat field correction (PFF) [24], uncorrected (UC); filtering: non-local means (NLM) [25,26], uncorrected (UC); contrast adjustment: gamma (G) [27], enhance contrast (EC) [28], uncorrected (UC); segmentation: multilevel thresholding: k-means clustering (cluster number) (KM(X)) [29,30], binary thresholding: default iterative intermeans (D) [31], Huang's fuzzy 2 (H) [32] (implemented by Schindelin), intermodes (IM) [33], isodata (ID) [31], Li's minimum cross entropy (L) [34][35][36], maximum entropy (MaE) [37], mean (M) [38], minimum error (MiE) [39], minimum (Mi) [33], moments (Mo) [40], Otsu's threshold clustering (O) [41], percentile (P) [42], Renyi's entropy (R) [37], triangle (T) [43], Yen's thresholding (Y) [36,44], ignore black voxel: iB, ignore white voxel: iW [45]. The source code of the binary thresholdings can be found in [45]. The different processing methods were written in the applied order from the left to the right and separated by "/". The combination of two binary thresholdings were indicated by "//". The image was copied twice before the segmentation process. One binary thresholding method was applied to each copy and both obtained thresholds were used to segmentate the original image. For example, the description Ref/RB/UC/G3.7/EC18/D-iBiW//R-iB means that the pathway bleach corrected by a reference image, lighting corrected with a low degree polynomial function, no filtering, gamma function with the input value 3.7, contrast en-Pharmaceutics 2021, 13, 1488 4 of 24 hanced with 18 % saturated voxels and thresholds determined with Default's threshold algorithm ignoring black and white voxels and Renyi's entropy algorithm ignoring black voxels. The details of the entire image processing evaluation and the different methods can be found in the previous publications [9,10]. Table 1 depicts the different matrices of batch T25E75 which were used for the in silico dissolution calculation. The corresponding matrices of the batches T50E50 and T75E25 can be found in the Appendix A (Tables A1 and A2). A sensitivity analysis was performed on the the third tablet of each batch. Bollmann and Kleinebudde (2021a,b) evaluated the image processing pathways for four binary mixtures (2 APIs, 2 excipients) in three different ratios each (12 batches). All pathways were analysed for different levels: API specific (all tablet images of batches containing the same API), combination specific (all tablet images of batches containing the same API and excipient), batch specific (all tablet images of one batch-same API, excipient and their ratio), and tablet specific. The pathways were ranked according to their desirability. The desirability was determined for the recovery rates of API mass, excipient mass and total mass and their respective ratios. The tablet desirability was calculated as geometric mean value of these six desirability values, each characterising a parameter of the virtual matrix. The API, combination or batch desirability, respectively, was calculated as the geometric mean of the geometric mean tablet desirabilities and a corresponding precision factor. A desirability of 1 is the optimum of the specified parameters. For more details on desirability assessment and image processing evaluation, see Bollmann and Kleinebudde (2021a,b) [9,10]. The final processed images of the best two pathways of the API specific (API1 and API2) and combination specific image processing (C1 and C2) were applied, as well as ones of the four highest performing batch specific image processings (B1-B4). Additionally, pathways with tablet specific desirabilities from 0.1 to 1 (increment 0.1) were used to analyse the sensitivity of the desirability (T1-T11).

Matrices of the In Silico Dissolution Calculation
Images of two batch specific pathways of the first and the fifth tablet were used to investigate the precision of the calculation of images of one batch. Furthermore, virtual tablets were designed by the TD-module and produced by the PAC-module to compare the dissolution profiles of the virtual tablets with the processed images of real tablets. The virtual tablets were produced in two different ways. API and excipient were either distributed randomly within the matrix (D) or random seeds were distributed and grown until the target tablet mass was obtained (G). Four different virtual tablets were designed for batch T25E75, both two with distributed API and excipient (D, D2) and two with seeded and grown compounds (G, G2). Table 2 depicts the input parameters of the dissolution calculation. The voxel edge length was 28 µm like in the images. The number of voxels in each dimension was set to the actual image dimensions. The dissolution profiles were predicted for 12 h with samplings every ten minutes. The mixing parameter was used to take the convection of the paddle stirrer into account. The identification numbers of the compounds were set to the actual grey values of the compounds in the image stack. The density of the phosphate buffer was determined by weighting a tarred volumetric flask. The molar mass of the compounds was taken from the certificates of analysis. The contact angle of ethyl cellulose was determined using a drop shape analyser (DSA 100, Krüss GmbH, Hamburg, Germany). The pore size was estimated after consultation with the developer of the F-CAD software package. The viscosity [46] and tension [47] of the phosphate buffer were taken from the literature.

Parameters of the In Silico Dissolution Calculation
The tablet porosity was calculated with Equation (1) with φ tablet as the tablet porosity, m total as total mass of the tablet, D as diameter of the tablet, h as tablet height, m API as API mass, ρ API as density of the API, m ex as excipient mass, and ρ ex as density of the excipient. It was used to determine the virtual porosity of the excipient for the dissolution simulation calculated by Equation (2) with φ ex as assumed excipient porosity for the dissolution calculation, S as number of slices of the image stack, φ tablet as tablet porosity, φ Image as porosity of the virtual tablet in the image, Vx size as edge length of a voxel, Vx API as number of voxels of the API, Vx ex as number of voxels of the excipient, and h as tablet height. The porosity of the virtual tablet in the images was calculated by the PAC-module.  The dissolution profiles were statistically analysed using python 3.8 [48]. The first parameter was the f 1 factor calculated using Equation (3) with f 1 as the difference factor of the calculated dissolution profiles, R t as cumulative percentage dissolved of the reference and, T t as cumulative percentage dissolved of the tested profile. The f 2 values were obtained by applying Equation (4) with f 2 as the similarity factor of the calculated dissolution profiles, R t as cumulative percentage dissolved of the reference and, T t as cumulative percentage dissolved of the tested profile. The calculated dissolution profiles were compared to the experimental determined dissolution profile [49]. For the calculation, either all data points of the first 12 h of the experimental dissolution testing or all data points up to the complete drug release (if < 12 h) and the corresponding ones of the predicted profiles were used.

Kinetic Analysis
The profiles were analysed regarding their kinetics by using the Higuchi and Korsmeyers-Peppas plot [15][16][17][18]. The linear regression line was calculated for all data points within 20% to 60% drug release. If less then 4 data points were in this range, the profile was excluded from the kinetic analysis. The equations of the linear regression lines and the adjusted R 2 of both plots were used to determine the kinetics. Additionally, slope and intercept of the regression line of the simulated dissolution profile were statistically tested against the regression line of the experimental dissolution profile using t-tests [50,51].

Determined Parameters
The results of the compound characterisation to determine the parameters for the dissolution simulation can be found in Table 2. The helium density of theophylline and ethyl cellulose was already determined for a previous publication [9]. The contact angle was determined as 93.2 ± 1.0°. The DS-module assumed that the tablets could not be wetted because the contact angle was above 90°. However, it was known from the experiments that the tablets were wetted. The DS-module only used the contact angle of the individual compounds, but not the one of the mixture. This may led to the wrong assumption. Therefore, the contact angle of ethyl cellulose was set to 89°which was the closest value to the determined parameter and led to wettability in the prediction.

Computation Performance
Although the dissolution testing lasted up to 72 h, the prediction of the dissolution profiles was only performed for the first 12 h due to the calculation time (approximately 13 h calculation time for 12 h dissolution time when running two simulations simultaneously, each using one GEFORCE RTX 2080Ti grafic card). The experimental dissolution testing is usually performed with six individual tablets. Therefore, the experiment was approximately three times faster than the in silico prediction of the DS-module. The calculation speed is dependant on the number of CUDA cores, the processor speed, and the speed of the memory. Therefore, it would be possible to reduce the calculation time with different hardware. In addition to the stand-alone version used in this study, there will be a cloud version available soon. This should accelerate the calculation as well. The deviations of the predicted profiles from the experimental profile after 12 h were estimated from the predicted profiles and the calculated parameters of the plots. Figure 1A depicts the experimental determined dissolution curve and different predicted dissolution profiles. There was a broad scattering of the predicted profiles using processed XµCT images as matrix with a tendency to higher deviations with a lower desirability of the matrix. The two dissolution profiles of the virtual produced tablets (T25E75-03_D and T25E75-03_ G) differed considerably from the experimental one and the XµCT image based predicted profiles. Figure 1B depicts that the virtual tablets produced by the PAC-module had a low performance. Both were considerably faster dissolving and reached a plateau at about 42% or 75% of released API after one or three hours, respectively. This indicates that the API was not percolating the virtual matrices. The experimental dissolution reached a plateau after approximately 72 h (data not shown). The API was released completely. Therefore, the API was percolating the tablet. The XµCT image based matrices seem to represent the percolating API structure. However, the kinetics of the dissolution profiles seem to be different from the experimental one. The graphs displayed in Figure 1C did not match the experimental profile but they were close to it. There was no clear increase in the performance visible for a higher specification of the image processing (API vs. combination vs. batch specific). This was surprising because a more specific image processing should increase the desirability. The sensitivity analysis Figure 1D depicts the broad scattering of the predicted profiles described above. Although the tendency was not clearly observed for the different specification levels, there was a dependency of desirability and prediction performance. A higher desirability led to a predicted profile closer to the experimental one. However, there were exceptions like T25E75-03_ T7. Therefore, there has to be a factor for the image quality which was not full described by the desirability. In spite of high performance of the API specific pathways, a combination specific image processing is still recommended because it increases the probability of higher desirabilities. Therefore, the recommendations of Bollmann and Kleinebudde [10] are still valid.  Table 3 displays the statistical analysis of the predicted profiles compared to the experimental one. Most f 2 -values of the XµCT image based predicted profiles were above 50 and most f 1 -values were lower than 15 which is the originally used limit for the similarity of dissolution profiles [49]. The most common used limit is a f 2 factor above 50 [12][13][14] which was obtained by even more virtual matrices. Although there were some exceptions (e.g., T7), the f 1 values were higher and the f 2 values lower for lower desirabilities. This confirmed the observed tendency. Table 3. The calculated results of the comparison of the predicted dissolution profiles and the experimental one. The first row in a cell is the value of the Higuchi plot, the second row represents the value of the Korsmeyers-Peppas plot. 'NA' indicates that a value was not calculated because either there were too few data points or it was not necessary, e.g., if the slopes were significantly different, the p-value of the axis intercept was not calculated.  Figure 2 illustrates the background of these exceptions. The f 1 and f 2 value of the virtual tablets T1 ( Figure 2B) and T2 ( Figure 2C) were similar. The pattern matched both well to the raw image ( Figure 2A). The slight differences of T1 and T2 were caused by a higher amount of assumed small particles (T2 > T1). The dissolution curve of T1 matched better to the beginning of the experimental curve, T2 was superior for the last part of the profile. Although the virtual tablet T7 had only a desirability of 0.506, the f 1 and f 2 values were superior. The pattern of the API was close to the one of T1. Therefore, the dissolution profile of the matrix matched the first part of the experimental curve. The recovery rates were 84.3% (API) and 99.9% (excipient) [10]. Some parts of the API particles were misallocated as excipient, some parts of the excipient particles were assumed as air. The overestimated pore size accelerated the calculated drug release. Therefore, most parts of the calculated profile matched the experimental curve. However, it is obvious, that a longer calculation time would increase the deviation because the difference of the profiles increased with the time. The cross-section of T4 explains the low performance of a matrix with a comparable high desirability ( Figure 2E, 0.801). The recovery rates were 110.5% and 98.3% [10]. The overestimation of API removed almost all pores and led to artificial fines at the edges of the matrix. The DS-module corrected the lower porosity, as described in Section 2.6. The artificial fines led to a faster drug release and a considerable difference to the experimental curve. This emphasises the crucial importance of a visual verification of the image processing. The desirability increases the probability of a suitable virtual matrix but does not guarantee it. It was only one matrix used for each desirability except for the desirability of 1. If all matrices would be used to calculate the dissolution profile, it would likely confirm the hypothesis, that most matrices with a higher desirability are superior. However, this needs to be investigated. This may provide information how the desirability calculation could be improved to reduce the user dependency of the visual verification. The profiles of the virtual tablets produced by the PAC-module were inferior. The f 1 and f 2 values of the virtual tablet using the distributed compounds were superior to the one using the grown compounds which was caused by the lower plateau. If the prediction was performed over a longer time period, the results would shift. Figure 3 depicts the differences of the virtual tablets created by the PAC-module and the image of the tablet. The XµCT based matrix ( Figure 3B) had a comparable distribution of API and excipient as the binned raw image ( Figure 3A). The virtual matrix produced with the PAC-module distributing the compounds ( Figure 3C) had a completely different appearance, while the virtual matrix obtained by seeding and growing ( Figure 3E) appeared similar. Both virtual matrices of the PAC-module had the overestimation of the pore size in common. Because both matrices reflected the full porosity of the tablet, there was no internal correction by the DS-module. Therefore, all API particles which were included by the excipient ethyl cellulose were not released. The fast increase in released drug and the plateau of the dissolution profiles were the product of too large assumed pores and a non-percolating API pattern of the virtual tablets. The particle size was designed by clustered voxels but it was not possible to check the particle size distribution and number of particles within the virtual tablets in the PAC-module. The seeded and grown virtual tablet was built by visual perception. The design of the seeded and grown virtual tablets was user depended because the user decided the seeded mass of each compound (determines the number of particles) and the number of growing steps (determines the particle size). Therefore, the designed virtual tablets may were not the optimum. Artificial virtual tablets produced by the PAC-module, which kept the ratio of API and excipient constant but reduced the interparticular porosity to 0 ( Figure 3D: distributing compounds, Figure 3F: seeding and growing compounds), changed the calculated dissolution profiles. Figure 1B depicts the differences. The internal correction of the porosity (the difference of the determined porosity and the interparticular porosity of the virtual tablet is assumed as intraparticular porosity) was able to mimic a percolating API cluster and the drug release was slowed down. Nevertheless, both profiles were inferior to the XµCT based virtual matrices and the calculated drug release was considerably faster than the experimental one. The matrix with the distributed compounds was superior to the seeded and grown one. This was caused by the small particle size and the higher number of particles.

Matrix
The Higuchi plot of the experimental dissolution profile indicates, as expected, a square root of time kinetic. The Korsmeyers-Peppas plot confirm the assumption as the parameter n was 0.48 which is close to the 0.45 proposed for a cylinder and Fickian diffusion as the transport mechanism of the API [12]. The slopes of both all XµCT image based matrices and the virtual tablets produced by the PAC-module were significantly different from the experimental determined drug release profile. The XµCT image based matrices obtained values for the parameter n in the range of 0.60 to 0.75 in the Korsmeyers-Peppas plot, while the parameter n was always lower than 0.4 for the virtually produced tablets by the PAC-module. The DS-module predicted an anomalous transport of the API for XµCT based virtual matrices [12]. Although the similarity of many XµCT image based matrices was given in accordance to the guidance of the FDA and EMEA [13,14], the profile itself differs significantly. Costa and Lobo pointed out some weaknesses of the f 1 and f 2 factors which were confirmed by the results presented [12]. Both kinetics and similarity did not match for the virtual tablets of the PAC-module. It is not recommendable to use the similarity of dissolution profiles to validate in silico predicted dissolution profiles. The kinetics and mechanisms of the drug release should be analysed as well. Figure 4 depicts the mean dissolution profiles of three tablets. The standard deviations of the experimental mean dissolution profile were low. The two mean curves of the XµCT image based matrices obtained similar standard deviations. The standard deviation of the virtually produced tablets were slightly increased but acceptable as well. The shape of the profiles did not change. Therefore, the calculations of the DS-module were reproducible and reflected the experimental deviations of different tablets. The different shapes of the predicted dissolution profiles seem to be systematic and should be investigated in more detail to improve the software package. Table A3 of the Appendix A depicts the data of the statistical analysis of the profiles. The results were in line with the results described for tablet 3. All slopes were highly significant different from the experimental curve.
XµCT based virtual matrices and the calculated drug release was considerably faster than the experimental one. The matrix with the distributed compounds was superior to the seeded and grown one. This was caused by the small particle size and the higher number of particles.
The Higuchi plot of the experimental dissolution profile indicates, as expected, a square root of time kinetic. The Korsmeyers-Peppas plot confirm the assumption as the parameter n was 0.48 which is close to the 0.45 proposed for a cylinder and Fickian diffusion as the transport mechanism of the API [12]. The slopes of both all XµCT image based matrices and the virtual tablets produced by the PAC-module were significantly different from the experimental determined drug release profile. The XµCT image based matrices obtained values for the parameter n in the range of 0.60 to 0.75 in the Korsmeyers-Peppas plot, while the parameter n was always lower than 0.4 for the virtually produced tablets by the PAC-module. The DS-module predicted an anomalous transport of the API for XµCT based virtual matrices [12]. Although the similarity of many XµCT image based matrices was given in accordance to the guidance of the FDA and EMEA [13,14], the profile itself differs significantly. Costa and Lobo pointed out some weaknesses of the f 1 and f 2 factors which were confirmed by the results presented [12]. Both kinetics and similarity did not match for the virtual tablets of the PAC-module. It is not recommendable to use the similarity of dissolution profiles to validate in silico predicted dissolution profiles. The kinetics and mechanisms of the drug release should be analysed as well. Figure 4 depicts the mean dissolution profiles of three tablets. The standard deviations of the experimental mean dissolution profile were low. The two mean curves of the XµCT image based matrices obtained similar standard deviations. The standard deviation of the virtually produced tablets were slightly increased but acceptable as well. The shape of the profiles did not change. Therefore, the calculations of the DS-module were reproducible and reflected the experimental deviations of different tablets. The different shapes of the predicted dissolution profiles seem to be systematic and should be investigated in more detail to improve the software package. Table A3 of the Appendix A depicts the data of the statistical analysis of the profiles. The results were in line with the results described for tablet 3. All slopes were highly significant different from the experimental curve.   Figure 5A illustrates the predicted and experimental determined dissolution profiles of the third tablet of batch T50E50. The experimental dissolution profile reflects a slowly eroding ethyl cellulose matrix. The profiles of the XµCT image based matrices were much less scattering compared to the ones described in Section 3.3, but the drug release was in all cases considerably faster than determined by experiment. The decreased scattering reduced the importance of the desirability and the visual verification since the sensitivity of the DS-module for those differences was lower. The virtual tablets of the PAC-module had a percolating API cluster but the drug release was completed within 40 min in both cases ( Figure 5B). This was probably caused by the assumption of artificial large pores, as shown in Figure 3. There was no clear difference between the matrices based on the API, combination or batch specific image processing ( Figure 5C). Although all curves showed a considerably slower dissolution than the virtual tablets produced by PAC-module, the predicted dissolution profiles were considerably faster than the experimental one. There seem to be a slight tendency, but the sensitivity analysis indicates only a negligible impact of the desirability on the prediction ( Figure 5D). This supports the conclusion of Section 3.3. The desirability represented not all important parameters and the combination specific image processing is still recommended as compromise between effort and outcome.

Batch T50E50
The results of the statistical analysis are presented in Table A4 of the Appendix A. The f 1 and f 2 factors confirmed a dissimilarity for all predicted profiles and the differences of the slopes of the Higuchi and Korsmeyers-Peppas plot were always highly significant. The parameter n of the Korsmeyers-Peppas plot was 0.47 for the experimental profile which is close to the 0.45 proposed for a cylinder and Fickian diffusion as the transport mechanism of the API as well [12]. Therefore, the erosion of the matrix was slower than the drug release. The parameter n was between 0.7 and 0.8 for all virtual matrices based on XµCT images. The DS-module predicted again an anomalous transport of the API for XµCT based virtual matrices [12] but the results were closer to the zero order kinetic than to the Fickian diffusion. This indicates that the DS-module overestimated the erosion of the matrices. It was not possible to calculate the kinetics of the virtual tablets produced by PAC-module due to the fast predicted drug release. The performance of the DS-module was worse compared to Section 3.3. The systematic error needs to be investigated to get reliable predictions of dissolution profiles. Figure 6 illustrates the mean dissolution profiles of three tablets. As described in Section 3.3, the shapes of the profiles remained identically with small standard deviations in all cases. Table A3 of the Appendix A depicts the data of the statistical analysis of the profiles. The results supported the results described for tablet 3. All profiles were dissimilar and the slopes were highly significant different from the experimental curve. This supports the assumption of a systematic calculation error within the software package which has to be investigated and corrected.    Figure 7A displays the predicted and experimental determined dissolution profiles of tablet 3 of batch T75E25. The different virtual matrices based on XµCT images resulted in only small differences of the predicted dissolution profiles. The profiles matched the experimental curve well. The profiles of the virtual produced matrices by the PAC-module were considerably faster releasing the drug and different from the other profiles. Figure 7B illustrates the considerable difference of the virtual tablets designed by the PAC-module. Figure 7C depicts that there was no relevant difference between virtual matrices based on different specific XµCT image processing pathways. Figure 7D displays the sensitivity analysis. There was no clear pattern visible. This underlines the conclusions drawn in Section 3.4. The desirability described not all relevant parameters. The recommendation to use the combination specific pathways stands because the API specific image processing led to huge variabilities within the image processing.  Table A5 of the Appendix A depicts the results of the statistical analysis of the profiles. The virtual produced tablets by the PAC-module were dissimilar to the experimental curve and a kinetic analysis was not possible due to the fast drug release. This was probably caused by the assumption of artificial large pores, as shown in Figure 3. All virtual matrices based on XµCT images were similar to the experimental profile. The f 1 -value was always lower than 9.5 and the f 2 -value always above 58. They all passed the originally determined limits ( f 1 < 15 and f 2 > 50) [49] and fulfil the requirements of the FDA and EMEA [12][13][14]. There seem to be no dependency of the similarity and the desirability left because there was no trend visible. This would reduce the effort of the image processing and the visual verification considerably. The parameter n of the Korsmeyer-Peppas plot was 0.74 for the experimental curve and indicates anomalous transport of the API [12]. The matrix was eroding. Therefore, the kinetic depended on the erosion and Fickian diffusion. The parameter n of all predicted matrices were statistically significant different from the experimental curve but the p-value increased considerably compared to Sections 3.3 and 3.4 (from range 10 −59 to 10 −113 to range 10 −4 to 10 −13 ). The slopes of the Higuchi plot mostly were not statistically significant different from the one of the experimental profile but the intercepts were always statistically significant different. Figure 8 illustrates the mean dissolution profiles of three tablets. The shapes of the profiles remained the same by low standard deviations. Therefore, the results were reproducible. The results of the statistical analysis (Table A3 in  virtual matrices based on XµCT images were similar to the experimental profile. The f 1 -value was always lower than 9.5 and the f 2 -value always above 58. They all passed the originally determined limits ( f 1 < 15 and f 2 > 50) [49] and fulfil the requirements of the FDA and EMEA [12][13][14]. There seem to be no dependency of the similarity and the desirability left because there was no trend visible. This would reduce the effort of the image processing and the visual verification considerably. The parameter n of the Korsmeyer-Peppas plot was 0.74 for the experimental curve and indicates anomalous transport of the API [12]. The matrix was eroding. Therefore, the kinetic depended on the erosion and Fickian diffusion. The parameter n of all predicted matrices were statistically significant different from the experimental curve but the p-value increased considerably compared to Sections 3.3 and 3.4 (from range 10 −59 to 10 −113 to range 10 −4 to 10 −13 ). The slopes of the Higuchi plot mostly were not statistically significant different from the one of the experimental profile but the intercepts were always statistically significant different. Figure 8 illustrates the mean dissolution profiles of three tablets. The shapes of the profiles remained the same by low standard deviations. Therefore, the results were reproducible. The results of the statistical analysis (Table A3 in the Appendix A) were in line with the results presented above. The DS-module was able to predict the dissolution profiles of batch T75E25. In contrast to the results of Sections 3.3 and 3.4 the statistical significant deviations had no practical relevance if virtual matrices based on XµCT images were used. The sensitivity of the DS-module regarding the virtual matrices based on XµCT images decreased with increasing API content of the tablets. This indicates more effort and a critical image processing, as well as visual verification for lower API doses. The reason for this phenomenon may be the reduced impact of the excipient on the dissolution profile. The drug release of the fast eroding tablets (batch T75E25) was predicted well by the DSmodule. This indicates that the dissolution of the API was predicted correctly. However, the erosion of the matrix was probably overestimated for batch T50E50. This resulted in a too fast predicted drug release. If the Fickian diffusion through a matrix was the drug releasing mechanism, the scattering of the drug releasing profiles of different virtual matrices increased. Therefore, the image processing seems to be more crucial for matrix based systems. The significant different kinetics of the predicted curves and the good prediction of batch T75E25 suggest a systematic error of the Fickian diffusion calculation The sensitivity of the DS-module regarding the virtual matrices based on XµCT images decreased with increasing API content of the tablets. This indicates more effort and a critical image processing, as well as visual verification for lower API doses. The reason for this phenomenon may be the reduced impact of the excipient on the dissolution profile. The drug release of the fast eroding tablets (batch T75E25) was predicted well by the DSmodule. This indicates that the dissolution of the API was predicted correctly. However, the erosion of the matrix was probably overestimated for batch T50E50. This resulted in a too fast predicted drug release. If the Fickian diffusion through a matrix was the drug releasing mechanism, the scattering of the drug releasing profiles of different virtual matrices increased. Therefore, the image processing seems to be more crucial for matrix based systems. The significant different kinetics of the predicted curves and the good prediction of batch T75E25 suggest a systematic error of the Fickian diffusion calculation of the API. The predicted transport mechanisms were different (Fickian diffusion vs. anomalous transport). This needs to be investigated and fixed.

Batch T75E25
The virtual tablets produced by the PAC-module behaved considerably different from the virtual matrices based on XµCT images and always were inferior. The virtual tablets of the PAC-module did not reflect the structure of the real tablets and led to considerably wrong profiles. Therefore, the design of those virtual tablets has to be revised critically and has to be changed until real tablets can be represented correctly. It may be possible to build similar tablets with the used PAC-module version with high effort but the benefits of the in silico tool would be lost. Although the software performed well for the eroding batch, it has to be improved to be a useful tool for the dissolution prediction in general.
With the current software version, the prediction performance is expected to be high for virtual matrices based on XµCT images of tablets which dissolution is determined by the solubility of the API (fast eroding systems, immediately release, low solubility of the API or large crystals). The data were reported to CINCAP so that future versions of F-CAD might be able to predict other systems, such as (slow eroding) matrices correctly with virtual matrices based on both XµCT images and the software itself.

Conclusions
The virtual tablets produced by the the PAC-module of F-CAD did not represent the real tablets. The predicted dissolution profiles were significantly different from the experimental profiles and the XµCT image based predicted profiles. It is mandatory to critically revise the design of virtual tablets within the software. In spite of a sufficient prediction of the dissolution of eroding matrices (batch T75E25) using XµCT image based matrices, the DS-module of software was not able to predict the dissolution of tablets with Fickian diffusion as drug releasing mechanism (batches T25E75 and T50E50) correctly. A systematic error in the calculation of these systems was observed which needs to be investigated. It is not recommended to judge the performance of in silico tools predicting dissolution profiles only by the f 1 and f 2 similarity [49] but also by the release mechanism and the respective kinetics to avoid misinterpretations (depicted with batch T25E75). The presented work emphasises the need of external validation methods of in silico tools. If an in silico tool designs virtual dosage forms and the software parameters are fitted until the predicted curve matches the experimental one, an optimisation within a certain design space is possible but no complete in silico development will be achievable. Additionally, the software may seem to be more accurate than it really is because some potential prediction errors could be concealed by over-fitting. The parameters have all to be determinable by the individual compounds and the in silico tools have to calculate the resulting physicalchemical properties of the mixture correctly. Otherwise, a complete in silico formulation development will not be possible. This is still a long way to go. Virtual matrices based on XµCT images provide the possibility to represent the real structure and distribution of pharmaceutical dosage forms. They can be used to validate both the virtually produced dosage forms of in silico tools and their performance. They provide the possibility to identify systematic errors and to optimise the performance of in silico tools. If the predicted dissolution curve considerably deviates from the experimental profile of the imaged dosage form, the known calculations of the tool may serve as starting point to investigate the differences. Furthermore, the use of virtual matrices based on XµCT images in in silico tools may help to increase the understanding of different mechanisms and interactions. XµCT image based matrices were successfully used to rate the performance of an in silico tool. The validation methodology presented was feasible to distinguish between a reliable and poor performance and to detect prediction errors of an in silico tool. The whole concept and work flow from the imaging and image processing [9,10] to the upload and calculation of virtual matrices based on XµCT images with an in silico tool worked well. However, there are still some improvements mandatory to use the presented methodology as a validation method. The sensitivity of the tested software package F-CAD regarding the desirability of the processed XµCT images decreased with an increasing API content. Therefore, both the image processing and visual verification becomes more critical with lower drug loads. One possible reason for this phenomenon could be the calculation issue of the software. Another important reason could be the pattern of the API and excipient within the virtual tablet. Although the overall pattern was similar, there were considerable differences regarding the particle number and size. The pattern was rated by visual perception and not included in the calculation of the desirability. However, the pattern seems to be an important parameter because the performance of the virtual matrices was not completely linked to the desirability. Therefore, either the visual verification has to be done carefully for all matrices with a high desirability (leads to a strong user dependency) or the pattern has somehow to be included in the desirability calculation. More investigations of this phenomenon are mandatory to understand the issue and to optimise the methodology.
Author Contributions: S.B. and P.K. contributed for the conceptualisation, writing-review and editing and methodology. S.B. contributed for the investigation, formal analysis, data curation, visualisation, software, validation, and writing-original draft preparation. P.K. contributed for the supervision. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The virtual matrices and the XµCT images of the nine presented tablets inclusive raw data of the XµCT measurement will be provided in a repository due to storage capacity issues. The link to the repository is available on request.

Acknowledgments:
The authors would like to thank Ashland and BASF for their support by providing the API or excipient, respectively. Furthermore, the authors thank CINCAP for providing the software F-CAD and Maxim Puchkov for the great software support. A special thank you goes to Thomas Bollmann, who supported the authors to translate their ideas to scripts to enable an efficient analysis of the data.

Conflicts of Interest:
The authors declare no conflict of interest.   NA  Table A5. The calculated results of the comparison of the predicted dissolution profiles and the experimental one. The first row in a cell is the value of the Higuchi plot, the second row represents the value of the Korsmeyers-Peppas plot. 'NA' indicates that a value was not calculated because either there were too few data points or it was not necessary, e.g., if the slopes were significantly different, the p-value of the axis intercept was not calculated.