Theoretical Modeling and Mechanism of Drug Release from Long

Article Cite This: Ind. Eng. Chem. Res. 2018, 57, 15329−15345

pubs.acs.org/IECR

Theoretical Modeling and Mechanism of Drug Release from LongActing Parenteral Implants by Microstructural Image Characterization Roberto Irizarry, Daniel Skomski,* Antong Chen, Ryan S. Teller, Seth Forster, Megan A. Mackey, and Li Li

Ind. Eng. Chem. Res. 2018.57:15329-15345. Downloaded from pubs.acs.org by DURHAM UNIV on 11/19/18. For personal use only.

Merck & Co., Inc., West Point, Pennsylvania 19486, United States

ABSTRACT: This work develops a mechanistic understanding and a simplified model to understand changes in drug release profiles of an active pharmaceutical ingredient in a long-acting parenteral polymer formulation from changes in its microstructure determined using X-ray computed tomography (XRCT) imaging. The system studied is an implant composed of a solid dispersion of crystalline drug microdomains embedded in a polymer matrix and fabricated by hot melt extrusion. We conduct material characterization of these devises using microCT and introduce, for the first time, imaging of such pharmaceutical devices by means of nanoCT, showing finer structures that cannot be captured by ordinary CT methods. We propose a simplified theoretical mechanistic model that estimates the drug release by extracting relevant microscopic features that modulate the release rate from microCT images to construct an idealized geometry. Using this analysis, we found that although the primary mechanism of resistance was originating from submicrometer features, which are below the resolution of the microCT, the main features that determine changes in diffusivity as a result of different process conditions were microscale pore volume and pore connectivity distributions. The parameters were extracted from the microCT images and in conjunction with the theory were able to explain changes in experimental drug release under different process conditions. These results demonstrate that the modeling strategy can explain the dominant release mechanism and estimate drug release of formulations, even when the equations are simplified due to many assumptions and a lack of subscale information. The proposed approach, due to its simplicity and reliance on the widely employed microCT technique, holds promise for guiding formulation development by providing a rapid initial assessment of batch consistency and process parameter selection before embarking upon lengthy and costly clinical trials.

1. INTRODUCTION In pharmaceutical formulation development there is strong interest in alternative delivery systems for the active pharmaceutical ingredient (API), that is, “the drug”.1,2 One particular area of interest lies with long-acting parenteral (LAP) formulations which achieve gradual drug release in patients over a period of many months or years.3 These formulations are intended to address the problem of poor patient adherence, which can be a major cause of therapeutic failure.3−5 Extended release formulations are intended to maintain clinically relevant doses of drug without requiring continued engagement on the behalf of the patient. Herein, an example of an LAP is studied which is a solid dispersion composed of crystalline API domains in a polymer matrix and © 2018 American Chemical Society

fabricated by a hot melt extrusion (HME) process. The aim is to devise a theoretical API release model based around the microdomain structure of the device extracted from X-ray computed tomography (XRCT) imaging. Historically,6,7 three release mechanisms had been proposed to explain extended release from matrix devices: variable diffusion coefficients, tortuosity, and pore constriction. It was concluded that pore constriction is a mechanism that was able to explain the very low diffusivities observed for this type of Received: Revised: Accepted: Published: 15329

June 22, 2018 October 9, 2018 October 17, 2018 October 17, 2018 DOI: 10.1021/acs.iecr.8b02806 Ind. Eng. Chem. Res. 2018, 57, 15329−15345

Article

Industrial & Engineering Chemistry Research

microdomain structures from imaging data in order to generate more realistic models which capture the entire interconnected device structure; called “network models”.20,30−35 Network models have the advantage that they are easy to implement and can tolerate some missing details during image analysis. Other methods consider the construction of irregular networks using training images to model the microstructure based on a crosscorrelation function.36 This approach has been used to study the impact of correlation length on transport regime using pore-network modeling for convection−diffusion systems.37 The analysis is extended to reactive transport.38 To guide formulation design in early pharmaceutical development it would be of enormous practical value to develop models that connect the microstructure to the drug release curves. Because microCT is used routinely in research laboratories, models based around this technique hold potential to accelerate the evaluation of devices fabricated with different processing conditions. The 3D imaging models provide a rapid initial read-out of projected release, which is much faster than the acquisition of experimental in vitro or in vivo release curves. Although 3D imaging has been used to model diffusion behaviors in a number of industries, to the best of our knowledge an image-driven modeling approach has not been applied to pharmaceutical long-acting implant formulations. We propose a set of image analysis techniques to extract key structural parameters from CT data and support our proposed model. We extract the drug domain as a porous network from 3D XRCT images and study the network connectivity to the external medium. The tortuosity of the network is studied statistically based on the geodesic distance from locations inside the drug domains to regions on the exterior of the implant. Similarly to that in ref 32 we conduct a skeleton extraction to set up the basis for understanding the geometrical and morphological properties, for example, the location and thickness of interparticle contacts, of the drug domains. Moreover, partitioning of the drug into individual components is performed using morphological operations. On the basis of the analysis, parameters for macro- and microlevel modeling are extracted. The microstructural details extracted from XRCT characterization are utilized with a simplified model that explains changes in drug release curves as a function of changes in the microdomain structure. Although CT has limited resolution, we demonstrate that the pore size and pore connectivity distributions can be used to generate insights into the dominant mechanisms that alter the release in formulations produced under different manufacturing conditions. Furthermore, these relative diffusivities can be calibrated with historical experimental release data to make the method semiquantitative and serve as a predictor of the expected release from formulations. Thereby, the method could offer estimations of release for new formulation batches for which no experimental release data are available. The paper is organized as follows: In section 2 we introduce the device fabrication and its characterization using XRCT. In section 3 we apply a macroscopic diffusional model to analyze experimental release curves of samples prepared by different process conditions. In section 4 we present a theoretical analysis and a minimalistic modeling scheme to predict effective diffusivity coefficients from the XRCT imaging data, thereby constructing a microscale diffusion model. The image analysis techniques developed for this model are also discussed

device. In order to understand and model drug transport mechanisms of controlled-release devices, theoretical models had been developed.8−16 Prior studies modeled release curves based on macroscopic diffusional and kinetic models. These models vary from each other in the diffusion mechanism that applies to the system under study (i.e., concentration- and porosity-dependent diffusion coefficients, dissolution-limited systems, etc.) as well as the geometry of the implant (i.e., spherical, cylindrical, etc.). The mathematical theories have been well established as means of understanding kinetic and thermodynamic release mechanisms and analyzing experimental release curves.12,13 One simplification to the diffusion equations is the Higuchi model in its different variations.14−16 This model has proven to be useful in many applications due to its simplicity and insights that can be gained from simple formulas.14 This type of model can also be misused when applied to systems for which the model assumptions and geometry do not apply.15 In the early stages of testing various modeling approaches, we implemented a macroscopic diffusion model to simulate release curves and to validate under which conditions the Higuchi assumption is valid for our system. The model tracks the moving boundary along the dissolution interface to model the transport of drug through open channels in the device. We found that for our system most of the dynamic dissolution behavior is approximated well by the Higuchi model. Only after 80% of the drug had released were deviations from the Higuchi model observed. Although macroscopic models such as the Higuchi model have been important tools in understanding diffusional and kinetic mechanisms that explain the observed release data, they do not offer information on how microstructural variations affect changes in diffusivity. The latter point is crucial because different process conditions can alter the microstructure of the device, and thereby, the release curves, which is particularly important for long-acting devices that need to achieve a target release rate over an extended period of time. For the modeling of micro- and mesoporous materials a large field of research has blossomed outside the area of controlled pharmaceutical release. Some examples of application areas are contaminant transport in clay materials for nuclear waste control, capillary entrapment of CO2, and petroleum extraction from rocks, among many others.17−25 Initially, microporous material models were studied assuming a random heterogeneous media26 or artificial network models,22,27 both generated artificially by computer algorithms following some underlying statistics. With advances in 3D imaging and image analysis techniques, these transport models have been enhanced by using the actual microporous domain information obtained from 3D imaging, for example, generated from XRCT and FIB-SEM, and then solving the underlying transport equations in the complex microstructures.24,25,27−29 These types of strategies vary from network models,23−25 to pore-scale models based on simulating the diffusion process with random walk particles,25 Lattice Boltzmann modeling,28 etc. Accurate predictions of the modeling approaches highly depend on the system studied, the sampling technique, the imaging instrument resolution, and the image analysis methodology. In terms of image analysis algorithms, many studies have been conducted to analyze these 3D images and a plethora of such works can be found in refs 17−19. Significant efforts have been devoted to develop algorithms to extract 15330

DOI: 10.1021/acs.iecr.8b02806 Ind. Eng. Chem. Res. 2018, 57, 15329−15345

Article


features (domain size, shape, etc.). XRCT was performed on the Xradia Versa XRM-500 scanner. The samples were imaged after fabrication without additional modification. The imaging software utilized is VersaXRM 10.7.3679.13921. X-ray source settings were 80 kV and 87 μA. The images utilized for the modeling were acquired at a 20× magnification. Reconstruction was performed with XMReconstructor software revision 10.7.3679.13921. Data representation was performed with ORS Dragonfly Pro v.2.0.0.0 (Zeiss license). The pixel size of data generated by this instrument under the scanning conditions tested is 0.7 μm. Note that while the terms “XRCT” and “microCT” are sometimes used differently (XRCT referring to a higher resolution method than traditional CT, capable of achieving resolution down to about one micrometer, as was used here), we use these terms interchangeably to distinguish microscale CT from the “nanoCT” method also discussed in this work. In 3D representations the cylindrical shape of the implants is clearly resolved (3D volume rendering in Figure 1a). The implants are not free of structural defects, but their abundance is quite low (≪1%) (2D cross sections in Figure 1b,c). In high magnification images, a contrast difference is observed between the API microfeatures and the polymer matrix. Zoom-in images at varying drug loadings (Figure 1d,e) provide compelling evidence that the API microdomains are resolved as bright dots in the images, whereas the surrounding polymer matrix appears dark. This contrast is intuitive because the physical density of crystalline drug is higher than semicrystalline polymer. Utilizing image analysis methods (outlined in section 4.2) the polymer matrix can be digitally removed from the images, leaving only the API microdomain structure (Figure 1f). A three-dimensional representation of the API microdomain structure is shown in Figure 1g. API domain properties extracted from these data sets were used for the drug release modeling. In Figure 1h a smaller data set comprising 10 cross-sectional slices is shown, which allows visualization of the degree of interconnectivity between API domains. For the API/polymer studied, the diffusion rate of the API through the polymer was determined by experiment to be negligible. Thus, for the creation of a diffusion model we consider that the API is not soluble in the polymer and the primary mechanism of drug release is diffusion throughout open channels that are formed as API is slowly removed over time through the interface with the solvent medium as shown in Figure 2.

in this section. The advantage of this strategy is that it integrates data with microlevel spatial resolution from XRCT, which is one of the main goals of this work. The second advantage of the model is that it provides a simple explanation of what drives changes in diffusivity coefficients. In section 5, nanoCT is introduced to generate higher resolution images of the device, showing additional features that cannot be captured by traditional XRCT, and provides an explanation regarding the limitations of traditional models. The conclusions are in section 6.

2. DEVICE FABRICATION AND 3D IMAGE CHARACTERIZATION In the present study the system of interest is one of a matrix implantable device, composed of API microparticles that are embedded and uniformly distributed in a polymer. The initial API size distribution and the extruder setup and process conditions (mechanical energy, temperature, etc.) determine the final microstructure. XRCT is utilized as an imaging tool to ascertain the macro- and microstructure of the implants. XRCT is an excellent imaging tool to extract microstructural data because it is commonly employed in industrial materials applications, can obtain three-dimensional images of devices nondestructively, and has sufficient density contrast and resolution to resolve API microfeatures embedded in the polymer matrix. We characterized samples made by two different process types. Generally the API and milled polymer are preblended and hot melt extruded to form rods of approximately 2 mm in diameter. Ten samples were prepared using the two processes varying some of the formulation parameter or process conditions with the objective of generating different microstructures. Two samples were generated using the process P1 at two different loading conditions. Then eight different samples using the process P2 were analyzed which changed in loading, temperature, and shear rate. Process 1 and process 2 (P1 and P2) refer to two different hot melt extrusion manufacturing equipment setups. The sample labeling, process type, and process details are summarized in Table 1. X-ray computed tomography (XRCT) was utilized to assess global uniformity (porosity, defects, etc.) and microstructural Table 1. Sample Labeling, Process Type, and Process Conditionsa sample no.

process

DL w/w%

S1_P1 S2_P1 S1_P2 S2_P1 S3_P2 S4_P2 S5_P2 S6_P2 S7_P2 S8_P2

P1 P1 P2 P2 P2 P2 P2 P2 P2 P2

40 45 40 45 45 45 45 45 40 45

temp

100 100 130 130

shear rate

3. MACROSCOPIC DIFFUSIONAL MODEL Prior to developing the image-based microstructural model, we sought to understand the general behavior of the release process by means of a detailed, macroscopic diffusional model which is outlined in this section and discussed in detail in Appendix A. Importantly, macroscopic models, such as the analogous Higuchi model, do not capture the effect of the microstructures. In practice, experimental release curves are utilized to fit the diffusivity of the macroscopic model to match the experimental data, thereby providing a value that is considered an estimation of the effective diffusivity coefficient. In this work, the macroscopic model is used to determine the effective diffusivity coefficient from the experimental release data and also to confirm the range of applicability of the Higuchi model for our system.

high low high low

a Temperature is given in units of °C. The shear rate is described in relative terms rather than an absolute quantity and is thereby dimensionless. P1 and P2 refer to “Process 1” and “Process 2.” P1 is a single-step extrusion manufacture, wherein a preblend was directly fed through a small twin screw extruder. P2 is a two-step extrusion, wherein a preblend was fed through a twin-screw extruder, pelletized, and then fed through a single-screw extruder.

15331


Article


Figure 2. Cartoon schematic of the domain structure of the implant. The implant is composed of microdomains of the active pharmaceutical ingredient (API, red) embedded in a polymer matrix (blue). Because diffusion of the API through the polymer is negligible, the drug release is primarily mediated by the formation of channels (gray) as drug from the surface dissolves and diffuses out of the implant and water penetrates deeper into the core of the implant. API domains that have connectivity to the surface of the implant are accessible.

Experimental evaluation of release curves is a standard technique to evaluate the behavior of a controlled drug-release device. It consists in placing the implant in a suitable ambient fluid and measuring the percentage of API released, W, as a function of time, t. For sake of model validation and calibration, we utilized the release curves generated for seven samples studied in this work. A typical release curve is shown in Figure 3 for two samples fabricated by the same equipment

Figure 1. X-ray computed tomography (XRCT) images of implant formulations containing API microparticles in a polymer matrix. (a) 3D volume rendering of a section of an implant containing 40% DL by weight with dimensions 2.25 mm by 1.95 mm. (b, c) 2D cross sections of the implant in panel a sliced perpendicular and parallel to the long axis of the implant cylinder, respectively. (d, e) Zoom-in 2D cross sections of implants with 10% and 20% DL by weight, respectively, showing the API microdomain structures which appear as bright dots. (f, g) Different representations of the API microdomain structure in 650 μm cylindrical segments inside the center of the implants. Only the API domains are shown (green) and the polymer is digitally omitted. In panel f, a sliced 3D representation is shown for an implant with 40% DL by weight. In panel g, a complete 3D representation is shown of an implant with 10% DL by weight. (h) A slice of thickness 7 μm from the representation in panel g, showing API domain interconnectivity. Scale bars are not provided for panels a and f−h because they are 3D representations. The cylindrical volume in panel a has a diameter of 2 mm. The cylindrical volume in panels f and g and the thin slice in panel h have a diameter of 0.6 mm.

Figure 3. Typical release curves, W(t), for S4_P2 (low shear rate) and S3_P2 (high shear rate).

but with different shear rates. We can see that changes in microstructures caused by different process conditions can greatly alter the drug release kinetics, even with an identical composition. The model we devised is an alternative to the Higuchi model and is used to demonstrate the region of validity for the Higuchi assumptions as applied for this system. These curves 15332


Article


Figure 4. Macroscopic modeling of the device. The macroscopic model serves as the precursor to the development of the microimage-based model. (a) In the region marked as “open pores,” the solid API is already depleted. In this depleted region, soluble API is transported by diffusion throughout open pores filled with ambient fluid. The diffusion coefficient is an effective diffusion coefficient controlled by the microdomains and nanoscale features. The API dissolving front is modeled as a moving boundary with radius h(t) that changes with time as the dissolution process proceeds. This interface is assumed to be at the equilibrium concentration. Finally, the ambient fluid is assumed to be a perfect sink, and this boundary is modeled imposing a zero concentration. (b) XRCT image, which shows how the dissolving front progresses from the periphery to the core of the device in the actual implant (18% of drug released). This front is modeled as a moving interface in the macroscale model.

follow the Higuchi-type release, in which the percentage of API released is proportional to the square root of time: W (t ) = k t

(1)

The linearity constant, k, which we call the release coefficient, is utilized in this work to measure the release performance of a device. We confirmed that when using side-by-side simulations between the detailed model in Appendix A and the Higuchi model, the Higuchi approximation was able to predict most of the profile with the exception of late stages (after 80% of total API released). The macroscopic diffusional model implemented in this work is based on a diffusional equation in 2D considering the implant’s symmetric cylindrical geometry, as is shown in Figure 4a. This assumption is confirmed using image analysis as shown in Figure 4b. The model domain is the diffusion of API through the open pores that are assumed to be filled with ambient liquid. This dissolving API front is modeled as a moving interface at the equilibrium concentration Cs. The external boundary is approximated as the sink boundary condition with a zero concentration. The details of this model are explained in Appendix A. The experiment’s effective diffusivity coefficient, Deff, is found by minimizing the difference between the model prediction of the release curve and the experimental data over the release time range: Deff = arg min|W (t ) − W exp(t )| D

Figure 5. Release coefficient vs effective diffusivity coefficient obtained from the macroscopic model: k vs Deff (k has units of mg ; Deff is dimensionless where D is the bulk diffusivity). day

Table 2. Release Coefficient and the Corresponding Effective Diffusivity Coefficients Obtained from the Macroscopic Model sample

release coefficient, k, (mg/sqrt(day))

S1_P1 S2_P1 S1_P2 S2_P2 S5_P2 S6_P2

5.39 8.98 1.97 3.72 6.38 2.96

experimental effective diffusivity coefficient, (Deff/D), dimensionless 2.4 8.0 4.0 1.5 3.5 9.0

× × × × × ×

10−03 10−03 10−04 10−03 10−03 10−04

4. MINIMALISTIC MICROSCALE THEORETICAL MODEL IN CONJUNCTION WITH IMAGE ANALYSIS TO PREDICT CHANGES IN EFFECTIVE DIFFUSIVITY COEFFICIENT In this section we introduce a microscale model that is sensitive to changes in device structure with the aim of explaining the cause of variations in the effective diffusivity coefficients from different processing conditions. The micro-

(2)

Plotting the effective diffusivity coefficient with the release coefficient it is clear that k ∝ Deff , (see Figure 5). A summary of the effective diffusion coefficients, and the corresponding release coefficients, for each sample is given in Table 2. 15333


Article


Figure 6. Idealization of the microstructure into tracks with cavities and necks in series (i.e., thin interparticle contacts that form nanochannels between primary particle domains).

Figure 7. Repeating unit for homogenization analysis highlighting primary domain dimensions that can be extracted from image analysis. The idealized geometry, shown here in an overly simplified schematic, is a series of connected cylindrical shapes.

approach must be implemented to address limitations of the method’s resolution. The minimalistic modeling approach proposed in this work has two objectives. The first objective is to understand if the limited features extracted from microCT images are enough to gain insight into the dominant structural features that modify the drug release when comparing different batches of formulations. The second objective is to determine if the model can be utilized to estimate changes in diffusivity as a result of alterations to the microdomain features imposed by different process conditions. We demonstrate in this section that the minimalistic model can satisfy both objectives. For the second objective we demonstrated that the model can be calibrated with data on macroscopic release curves to account for the unresolved features, for example, thin particle contacts, that are below the microCT resolution. The calibrated model extracts a semiquantitative prediction of the diffusion rate, Deff, of new formulations by means of 3D image analysis. This tool enables fast exploration of the design space without needing to generate a large number of release curves which have a very long lead time. The authors acknowledge that this approach

scale model, as the name suggests, differs from the macroscopic model in that it is sensitive to microstructural features that change between the different formulations tested. The microscale model discussed in this section is referred to as a minimalistic model because it incorporates certain assumptions and simplifications about the device structure. This simplified approach reduces the accuracy of the method, but results in a compact model. The model is able to provide an understanding of major release rate changes from a few salient features of the microCT image analysis. We also introduce the image analysis algorithms needed to support the model assumptions and to extract key features such as pore size distribution, pore connectivity, and coordination number which are used as model inputs to predict drug release. However, the microCT image does not have sufficient resolution to capture the thin contacts, or “necks,” between particles which are a main resistance for diffusion in our system. This resolution limitation hinders a direct application of microCT imaging data with first principle modeling to predict drug release (as discussed in section 5). Therefore, an 15334


Article

Industrial & Engineering Chemistry Research does not replace the need for 3D nanoimaging methods. Indeed, such advanced tools are discussed in section 5. Rather, the current approach is useful primarily due to its ability to provide a rapid read-out of projected release by making salient estimations and is reliant on an instrumental method that is commonly employed in industrial research laboratories. 4.1. Minimalistic Model Construction. A microscopic model for the effective diffusivity coefficient is proposed based on an idealized microstructure. The actual microstructure is idealized as radially directed and connected tracks. The tracks are composed of “cavities” and “necks” in series. The cavities represent particle agglomerates and the necks represent interparticle contact points between the aggregates (see Figure 6). The presence of thin particle contacts, or necking, is confirmed by the nanoCT data as discussed in detail in section 5. The authors acknowledge that there may be additional factors which could also contribute to the slow release such as the osmotic rupture of thin polymer layers between neighboring API particles. Interparticle constrictions are applied in this model because they present one effective mechanism to understand and model the slow release kinetics of this system. For the model construction it is further assumed that this microstructure is a repeating unit from which effective properties can be estimated using homogenization analysis39,40 (see Figure 7). The ansatz of homogenization theory provides general solutions for an arbitrary repeating geometry in the form of a complex set of partial differential equations. The justification for homogenization is that this work is intended toward the study of structurally homogeneous formulations. Indeed, a primary objective pursued in the pharmaceutical implant development of solid dispersion matrix devices is to reduce or eliminate macro-scale defects in the implants, producing structurally uniform formulations. Indeed, in the samples tested here differences in connectivity/ tortuosity between different regions in the same implants studied were found to be negligible. A small volume fraction of defects is present in some of the formulations, but did not impair the model’s ability to interrogate major pharmaceutically relevant trends between formulations. The authors acknowledge that the approach pursued is not valid for formulations exhibiting profound heterogeneities or for analysis of subtle changes in implant submicrometer structure. Rather, the goal is to create a compact calculation that can capture a sufficient amount of structural information to model significant variations in the drug release rate of the kind that would be meaningful to pharmaceutical development. This idealized microstructure allows utilization of a simple form of an effective diffusivity coefficient from homogenization theory:39,40 ⊥ Deff =

(L + l)D1D2 LD2 + lD1

To adapt these formulas to our case, we assume that the “layered” material is in the radial direction and the material for type 1 is the pore and the material type 2 is the constriction between connecting pores. In our case we need to generate estimates for D1 and D2 for the pore and constriction layers, respectively. We estimate the pore diffusivity coefficient as a bulk diffusivity coefficient and the neck’s diffusivity coefficient is proportional to the area reduction of the neck relative to the pore radius. iry D1 ≡ Dpore = D , D2 ≡ Dneck = Djjj zzz kR{

(LD1 + lD2) (L + l )

(5)

where R is the API particle radius, and r is the radius of the neck. The incorporation of structural features from image analysis into the equations is described in section 4.2 and the approach is outlined in Figures 6 and 7. Also, since there is a distribution of pore sizes and neck sizes, it is assumed that this distribution can be placed radially and the effective diffusion coefficient we want to calculate, which is the average radial diffusion coefficient, can be obtained by integrating eq 3 over a distribution of pore-sizes and neck sizes: Deff =

∫ f (L , l , R , r )

(L + l)D1D2 dL dl dR dr LD2 + lD1

(6)

Some of the assumptions of this idealized geometry are justified with further image analysis, as will be discussed in the subsequent sections. Nevertheless, the idealized geometry is still an oversimplification of real geometry composed of a very complex microstructure. In practice these cylindrical radii are not extracted directly, but correlated with parameters that can be extracted from image analysis. For example, the radius of the larger cylinders is associated with the pore volume and the radius of the smaller cylinders is associated with a nominal neck radius as well as the number of contacts of a pore with other pores. This bridge connecting theory and extractable parameters from XRCT image analysis is discussed in section 4.3. Finally, another simplification is that the pore volumes detected are not cylindrical but have a random shape and all the connectivities are not in radial direction but in random directions. 4.2. Image Analysis for Testing the Minimalistic Model Microstructure Assumptions and Determining Important Parameters Used by the Model. In this section we examine three image analysis methods that were used for understanding the relevant microfeatures that determine the release differences between formulations. A couple of preliminary methods were explored for mechanistic understanding. First, we examined the application of percolation theory to assess the degree of API domain connectivity. The technique extracts the API domain as a porous network from 3D XRCT images and examines their connectivity to the external media. Second, we applied a geodesic pathway analysis to understand the microchannels formed throughout the device. Similarly to that in ref 32, this geodesic analysis is conducted by extracting the skeleton of microchannel paths to set up the basis for understanding the geometrical and morphological properties of the API domain, for example, track thickness and throat locations. The geodesic distance aids in understanding the impact of network tortuosity on drug release. The two aforementioned methods (percolation theory and geodesic distance) were mainly applied as stepping stones toward gaining a mechanistic understanding of the drug release

(3)

This formula represents the diffusion coefficient in the direction of layered periodic materials, where material type 1 has dimension L and diffusivity coefficient D1 and material type 2 has dimension l and diffusivity coefficient D2. The formula for the parallel diffusivity coefficient is given by P Deff =

2

(4) 15335


Article


Figure 8. Visualization of implants through volume rendering. The polymer is digitally omitted from the representation. The gray features are the API microdomains. The black plane is merely an artificial backdrop to aid in visualization. Indicated are (i) the total API microdomain volume, (ii) the API domains that form connected pathways which reach from the interior to the exterior of the device, referred to as “accessible” API domains, and (iii) API domains which do not have connected pathways to the exterior of the device, referred to as “inaccessible” API domains. An analysis of inaccessible API domains provides an illustration explaining why only a small portion of API material can release at low drug loadings. A scale bar is not provided because the images are 3D representations. For each representation a thin cylindrical volume with a diameter of 0.6 mm was probed.

Figure 9. Geodesic distance from API voxels on the interior of a volume to the exterior (units in micrometers). Example transverse slices are extracted from implant samples with 30% (left) and 20% (right) drug loading. Note that the yellow color on the right panel represents an infinite distance (meaning that the domains on the interior are mostly inaccessible from the exterior of the device).

from the device. The final image analysis method, and the cornerstone of the simplified microscale-model, involves extracting the pore-size and connectivity matrix using a morphological erosion/dilation approach, which partitions the API into individual subcomponents. On the basis of the analysis, parameters for macro- and microlevel modeling are extracted. Note that although various XRCT image analysis techniques were discussed in ref 19, we did not have the need to incorporate them in our approach. The image analysis was conducted using programs implemented in MATLAB with 40 consecutive slices from the cylindrical implant microCT data sets. 4.2.1. Testing Connectivity As a Function of Percent Drug Loading (Percolation Theory Image Analysis). A 3Dconnected component analysis is performed on the API clusters, which are interconnected to the external medium and form a network of pathways for dissolution. For image analysis of complex microstructures, image segmentation into the principal components (e.g., API and polymer) is a challenging problem and a topic of continued research in the image analysis community. In our case, we opted for segmentation between the background, the polymer, and the API by

performing Otsu’s optimal thresholding method.41 We implement primarily three different classes: polymer, API, and an artificial background class. As a result, as determined by the density of individual materials, voxels in the lowest intensity range are assigned to the background class and voxels in the highest intensity range are assigned to the API class, and the remaining voxels with medium intensity belong to the polymer class. This analysis was translatable to a variety of different samples and able to describe the important microstructural changes induced by process variations. Binary masks in 3D for the three classes are then saved to facilitate further analysis. In this analysis, the foreground mask, i.e., the combination of the masks for API and polymer, is dilated42 by one voxel, and then subtracted by the original foreground mask, which generates a “shell” mask. Therefore, all the API clusters that are interconnected to the shell are considered accessible from the outside. Correspondingly, the API clusters that are disconnected are considered isolated and inaccessible. Figure 8 demonstrates two samples (30% and 20% drug loading) for which the accessible and isolated API clusters are visualized via 3D volume renderings. It can be seen that a substantial portion of API clusters are inaccessible with 15336


Article

Industrial & Engineering Chemistry Research 20% DL. The lack of substantial domain connectivity is consistent with corresponding drug release data, which indicates miniscule drug release. These examples demonstrate the utility of microimaging toward gaining a qualitative understanding of the microstructural effects that guide release from the devices. While a range of drug loadings was investigated in the development of the model, one key assumption of the theoretical analysis is that the drug loading is sufficiently high for the domain connectivity in the implant to be above the percolation threshold. From a practical standpoint this assumption is reasonable for the pharmaceutical formulations considered in this work. In solid and uniform matrix-based implants with negligible intrapolymer diffusion, significant quantities of drug are trapped in the formulation unless the microparticle network is fully connected to the implant surface and the external medium. A primary aim of the formulation is to maximize the duration of the release, hence requiring all of the API in the formulation to become fully accessible toward the dissolution (since the drug diffusivity directly through the polymer matrix is miniscule). Thus, formulations with network structures above the percolation threshold are the focus of this work because they were experimentally observed to exhibit complete release of the API. 4.2.2. Testing the Radial Pathways (Geodesic Pathway Image Analysis). Having gained basic understanding of the release using a percolation threshold analysis, we utilized a geodesic model which identified and helped understand the radial pathways for diffusion. On the basis of the 3D connected components, we are able to compute the length of the shortest path from a voxel inside an API domain to the exterior medium through the connected API pathways, that is, the voxel’s geodesic distance43 to the exterior shell. Figure 9 shows a transverse slice from implant samples with the voxel’s geodesic distance shown in a colormap. The colormaps illustrate that when the API voxels are inaccessible, as in Figure 9 on the right, their geodesic distance is infinite. On the basis of the distance measures, for each voxel in the API domain, the tortuosity τ is calculated as τ = G/D

Figure 10. Examples of API pathways extracted from geodesic image analysis, where arbitrary voxels inside the API domains are sampled and shown in red color, and their shortest pathway toward the exterior of the implant is extracted. The pathways are shown via volume rendering in orange color, and the contact points of the pathways are shown through surface rendering in the form of dots in gold color.

parameters extracted for implant images from various compositions and conditions. It is found that the tortuosity is generally low (1.1−1.2) for all samples. Also, lengths between domain bodies and constrictions were estimated (L and I, respectively). 4.2.4. Extracting Pore-Size and Connectivity Matrix Using Erosion/Dilation Technique. Although the percolation connectivity and geodesic pathway analyses discussed in the prior sections are useful in confirming some of the theoretical assumptions, further analysis would be needed to extract detailed geometrical features, for example, a quantitative estimate of API particle radius R and neck radius r based upon partitioning of the API domain network. Since the API pathways in the implants with high drug load tend to be fully connected, the partitioning of the entire API network into individual units is not a trivial task. Our approach allows us to partition the API network, which is actually a sizable supernetwork of interconnected paths, into smaller disconnected domains where the size and shape of the domains can be studied. Herein we perform a series of morphological operations to the binary mask of API IM starting with

(7)

where G is the geodesic distance and D is the Euclidean distance between the voxel and the shell. With the accessible API pathways identified, we are able to extract the skeleton of the pathways, and calculate their properties, such as the path length. As the representation in Figure 10 illustrates, for any given accessible API cluster, the shortest pathway to the external medium can be visualized as a chain of connected API clusters, and the lengths measured of the pathways can be extracted (the perpendicular distance from the points on the skeleton to their closest track boundary on the pathway). Note that the examples are only showing the pathway with the shortest geodesic distance. In the actual implant, since the API clusters are interconnected in 3D in the form of a network, there are multiple pathways for an API cluster to reach the external media. The geodesic pathway analysis primarily serves to demonstrate how it is possible to extract qualitative information from microstructural imaging data to visualize the pathways of drug release from the device. 4.2.3. Extracting Quantitative Information from Image Analysis. In addition to describing the connectivity of domains and pathways for API diffusion, many additional key properties can be extracted through image analysis. Table 3 shows

IM̂ _binary = (IM _binary ⊖E)⊖E

(8)

where E is a binary 3 × 3 × 3 structure element and ⊖ represents the erosion operation. By performing the erosion operations, the necking regions of two voxels or less are removed, which disconnects the API network into pieces. For ̂ the eroded mask IM_binary , a 3D connected component analysis is performed, such that disconnected pieces are identified and ̂ labeled as IM_labeled . On the basis of the labeled mask, a series of morphological dilation operations are applied to recover the API domain approximately to its original size, where the reconstructed and labeled mask is obtained as IM̃ _labeled = (IM̂ _labeled ⊕ E) ⊕ E

(9)

where the element E is the same as the one used for the erosion, and ⊕ represents the dilation operation. Examples incorporating this analysis are shown in Table 4. 4.3. Application of the Minimalistic Image-Based Model Approach (Connection between the Idealized Geometry and XRCT Image Analysis). The approaches 15337


Article

Industrial & Engineering Chemistry Research Table 3. Properties of Implant Samples Extracted through Image Processinga nominal API mass ratio actual API load volume ratio isolated API volume ratio track radius mean Estimated tortuosity mean

S1_P1

S2_P1

S1_P2

S2_P2

S3_P2

S4_P2

S5_P2

S6_P2

S7_P2

S8_P2

0.40 0.41 0 1.70 1.15

0.45 0.43 0 1.63 1.14

0.40 0.35 0.002 1.76 1.28

0.45 0.39 0.001 1.77 1.19

0.45 0.42 0 1.40 1.13

0.45 0.39 0.001 1.70 1.21

0.45 0.40 0 1.66 1.16

0.45 0.38 0.001 1.73 1.22

0.40 0.34 0.002 1.70 1.30

0.45 0.36 0.002 1.73 1.27

a

Image processing as in sections 4.2.1 and 4.2.2. Note that the nominal % API loading is given on a by-weight basis, whereas the API quantities determined by imaging are calculated on a by-volume basis.

Table 4. Particle Size Statistics for Different Samples and Formulations sample

release coefficient, k, (mg/sqrt(day))

average pore volume, (voxels)

average coordination number/particle

S1_P1 S2_P1 S1_P2 S2_P2 S3_P2 S4_P2 S5_P2 S6_P2 S7_P2 S8_P2

5.39 8.95 1.97 3.72 5.60 2.81 6.38 2.93 1.54 3.71

440 364 946 842 222 553 372 587 619 701

3.51 3.64 3.21 3.34 3.22 3.25 3.45 3.21 2.95 3.10

Deff, p

i rnom 2nc ,p2 y (L + l)jjj V zzz k p { ∝D 2 i rnom nc ,p2 y Ljjj V zzz + l k p {

The values of L and I are assumed to be nominal values that can be extracted from the geodesic analysis (explained in section 4.2). The value of rnom is assumed to be an undetermined radius. This equation is further simplified by considering that for small constrictions, the second term in the denominator is much larger than the first term, which can be ignored as a first approximation. Therefore, eq 9 can be further simplified to Deff, p ∝

discussed in the aforementioned section allow extraction of qualitative and quantitative microdomain features from imaging data. However, application of eq 1 and eq 3 directly from the XRCT images is not possible because some of the parameters needed cannot be resolved within the resolution limitations of the images, as is discussed in section 5. For example the neck radius, r, is below the XRCT resolution (ca. one micrometer). This is a nontrivial limitation because a main resistance toward diffusion is the geometry of the constrictions. Furthermore, the pores have a complex structure that is not easy to describe with a unique cylinder radius, R, as was discussed in section 4.1. To solve these limitations, we follow a softer goal of approximating some of these model parameters using the pore size and connectivity data extracted using the erosion/dilation image analysis method. The goal is to generate a functionality that links the effective diffusivity with the pore size distribution and connectivity. Since there are many parameters that are difficult to describe, we call this simplified functionality the relative effective diffusivity, Drel eff. To develop a formula for this relative diffusivity, the following simplifying assumptions are made: First, we assume that the neck radius, r, in eq 5 can be approximated by a nominal neck radius, rnom, and coordination number for each particle, nc. We also assume that the nominal neck radius does not change significantly from sample to sample (at the current stage of the discussion, this is a model assumption that cannot be confirmed a priori). Second, we assume that the pore radius, R, is proportional to the pore volume (typically of irregular shape). This assumption may be partially justified by the geodesic analysis where an estimate of the pore length was consistent for most samples, so we can assume that the change in volume can be due to the lateral area. 2 rp ∝ rnom nc2, p

and

R p ∝ Vp

(11)

nc2, p Vp

(12)

This theory suggest that the release coefficient can be represented by the following general formula, Deff, p ∝

ncα, p V pβ

(13)

where α and β are empirical parameters determined from historical experimental release data, which serve to calibrate the model. Given this functionality for a single pore dimension and connectivity, to determine the effective diffusivity coefficient, following eq 6, this kernel is integrated over the particle numbers and coordination number. Deff ∝ Z XRCT(α , β)

(14)

where Z XRCT(α , β) =

Vmax

∫V ∫n

nc ,max

α jij nc , j zyz jj zz/Np jj β zzz V j=1 j k j { min

c ,min

i nα y f (V , nc)jjjj cβ zzzz dV dnc kV {

Np

≈

∑

(15)

If this proposed theoretical model describes the dominant release dynamics, the parameter ZXRCT is linearly proportional to the experimental diffusion coefficient (Deff ∝ ZXRCT) with a calibration constant to account for missing information like unknown neck radius. Since the effective diffusion coefficient is also related to the release coefficient by Deff ∝ k2 (see Figure 5), the linear correlation can be written in term of k2. We can further simplify the formulation using k instead of k2 (which is justified when the range of k is not too large and also is absorbed in the empirical parameters α and β). This is a more practical correlation since then the release curve can be predicted using eq 1 without a need to solve the macroscopic diffusion model. In this case we propose the following linear

(10)

These parameters are substituted in eq 1, and the effective diffusivity using the properties of a given particle, p, is 15338


Article

Industrial & Engineering Chemistry Research correlation to predict the release coefficient k with the parameters ZXRCT: k = aZ XRCT(α , β) + b

The results demonstrate that this simplified model can be utilized to understand variations in formulation performance as a function of changes in the microstructure. This theory helps in two ways. First, the model gives a very simple explanation of the involved mechanism. The theory proposes that larger pore sizes reduce the effective diffusivity due to a reduced pore-size/ neck-size ratio, while larger coordination numbers increase diffusivity due to an effective increase in pore-size/neck-size ratio. This simple understanding helps engineers in decision making during process design and optimization. Furthermore, when some experimental data is available like that shown in this work, a linear correlation like the one shown in Figures 11 and 12 can be utilized as a calibration curve between the experimental and the relative (i.e., model) diffusivity coefficients. This, in turn, allows computing anticipated release curves of new formulations without the availability of release data. We propose the work flow shown in Figure 13.

(16)

We plotted the equation and, indeed, k followed a linear trend with eq 12 as shown in Figures 11 and 12. As an initial test, in

Figure 11. Comparison of the model (eq 15) with experimental data for unoptimized values of α = 2 and β = 1 (k has units of mg ; ZXRCT day

has units of

1 pixel

). β

Figure 12. Comparison of the model (eq 15) with experimental data for optimized values of α = 4.0, β = 0.87.

Figure 13. Proposed work flow for the simplified microstructural model to determine anticipated drug release. The image analysis (section 4.2) involves extracting the pore-size and connectivity matrices using a morphological erosion/dilation approach, which partitions the fully interconnected drug network into individual subcomponents and then computes their connectivity to neighboring microfeatures.

Figure 11 we utilized the parameters approximated in eq 10 (α = 2, β = 1) to estimate the parameter ZXRCT(α,β) with a squared regression coefficient of 0.70. Optimized values for the empirical model exponents, α and β, in eqs 14 and 15 can be determined using the following simple procedure: a 2D grid of (αi,βj) is generated in a specified range. For each value of αi and βj in the grid a linear regression is utilized to generate the model corresponding to these values. The values that produce a model with the lowest squared error between the model and the data are selected as the optimal empirical parameters to complete the definition of ZXRCT in eq 15. In our case, instead of the grid analysis we implemented this procedure as a formal optimization.44 Figure 12 shows the correlation with optimized parameters (α = 4.0, β = 0.87), in which case the squared regression coefficient is 0.90. This correlation can be utilized as a calibration curve to predict release coefficients from XRCT analysis.

The minimalistic model helps explain observed changes from process conditions. First, we examine the difference in drug loading under identical process conditions. Figure 14 shows that the effect of going from 40% to 45% drug loading is mainly an increase in connectivity. The range of the microdomain size distribution does not shift significantly, but there is an increase in the number of interparticle contacts. Correspondingly, there is an increase in the effective diffusivity as shown in eq 10, which is confirmed by the experimental data. The models indicate a picture of release by which the mechanism of increased diffusivity from increases in the drug loading % for a given process is primarily due to a change in 15339


Article


Figure 14. Distribution of pore size and particle coordination number for two different drug loadings (40% and 45%) at the same process conditions (samples S1_P2 and S2_P2, respectively).

Figure 15. Distribution of pore sizes and particle coordination numbers with the same drug loading, but two different processes P1 and P2 (twinscrew vs twin-screw/single screw extrusion) (for samples S2_P1 and S2_P2, respectively).

number of contacts, rather than the increase in microdomain sizes. Second, we examine the impact of changing from a twinextrusion setup (P1) to a twin-extrusion/single-extrusion setup (P2) under identical drug loading. Figure 15 explains why process P2 has a much slower release than process P1 with the same drug loading. The reason is a decrease in the number of interdomain contacts and also an increase in domain size. These quantitative observations are corroborated by a qualitative visual picture of the formulations from XRCT imaging (Figure 16), which agree with the suggested mechanism (specifically, variations in domain size between samples). Originally, it was conjectured based on a qualitative analysis of the XRCT images that tortuosity might be the main difference between these processes, due to a circular orientation of domains in P2 which is absent in P1, but rather it is the larger number of drug microfeatures and the increased connectivity that determine the main changes in diffusivity (which require a quantitative readout of the 3D connectivity from the images to assess). Third, using the same drug loading and equipment setup (P2), we examined the effect of varying the shear force during the extrusion (Figure 17) and the extrusion temperature (Figure 18). Figure 17 demonstrates that transitioning from a low shear to high shear process, the increase in diffusivity is due to a reduction in domain size and increased number of contacts. Finally, Figure 18 shows that processing at 100 °C and at 130 °C produce similar structures in terms of domain size distribution and connectivity distribution. The model

Figure 16. XRCT image 2D cross sections of (a) process P1 and (b) process P2 (samples S1_P1 and S1_P2; compare to Figure 16).

15340


Article


Figure 17. Distribution of pore size and particle coordination number for two different shear rates for process P2 (samples S3_P2 and S4_P2 produced with high and low shear forces, respectively).

Figure 18. Distribution of pore size and particle coordination number for two temperatures for process P2 (samples S4_P2 and S6_P2).

Figure 19. Comparison between the model prediction using the work flow shown in Figure 13 and the experimental data. (a) Comparison between low and high shear at 100 °C and (b) comparison between low and high shear at 130 °C.

Although the main resistance for diffusion is below the XRCT resolution, it is shown in this section that for the polymer system studied, a minimalistic approach can explain differences in release between formulations. The method extracts pore coordination numbers from microimage analysis and uses assumptions to fill in missing details from the submicrometer structure. The microscale modeling approach serves as a quick guide during formulation design. It gives a clear insight into the microscopic structural features that modify release between the formulations, which can be summarized as slow diffusion inversely proportional to domain

predicts that for these conditions the diffusivities should be very similar, which is confirmed experimentally. Lastly, the results of the minimalistic microscale model can be calibrated with historical release data from the macroscale model to generate anticipated release curves for new formulation batches. Figure 19 shows typical results of the minimalistic model obtained using the work flow described in Figure 13. As can be seen by these results, this simplified model using limited information from XRCT can predict release curves with less than 10% prediction discrepancies against experimental data. 15341


Article


application of pore-scale models generates artificially large effective diffusivity coefficients. As seen in this side-by side comparison, the pore volumes identified using image analysis from microCT are actually a composite of API particles and the spaces between them and the polymer. Another important feature of the nanoCT is that finer structures that can be associated with thin particle contacts/necks can be identified and measured as shown in in Figure 21. The images prove the presence of thin constrictions between particles, which were utilized for the modeling construction as a model simplification (outlined in section 4).

size and directly proportional to the number of domain contacts and contact radius.

5. INCREASING RESOLUTION WITH NANOCT AND ASSESSMENT OF OTHER MODEL STRATEGIES The minimalistic model developed in the previous section resulted in a simple and compact theory in which changes in the release coefficient, k, due to different process conditions can be estimated from the pore size and coordination number distribution determined by the erosion/dilation image analysis of XRCT images (see eq 15 and eq 16). These methods are advantageous because they allow assessing long-acting implants in a rapid manner and are reliant primarily on microCT imaging, which is a technique that is widely adopted in industrial and academic research laboratories for nondestructive 3D imaging. However, microCT cannot resolve submicrometer features which are important in controlling the rate of diffusion in matrix-based implants. The theoretical analysis developed in this work overcomes the resolution limitation through model calibration from historical release data, resulting in a simplified model to predict drug release from different process conditions that can be applied to new batches of formulations for which experimental release data are not available. Nevertheless, it would be desirable to improve analysis of the samples by incorporating 3D nanostructural data to enable a direct application of pore-scale models toward the release modeling. To understand submicrometer structure, we characterized an implant using nanoCT and conducted a direct comparison to microCT data. To the best of our knowledge, this is the first time that a pharmaceutical extended-release implant has been characterized using nanoCT. These images are shown in Figure 20, which demonstrate the massive 10-fold improvement in resolution of nanoCT compared to microCT. The resolution limitation in the microCT explains why a direct

Figure 21. Thin interdomain contacts (“necks”) form restricted channels between particles, as resolved using nanoCT. With the nanoCT we are able to explore in greater detail the neck structure between pores. For examples, the nanoCT is able to identify the necks shown with red arrows. The measured diameter of those necks is about 600 nm, below the resolution of the microCT which can only resolve microscale features.

The nanoCT data provides a key piece of evidence demonstrating that XRCT does not have the resolution to apply first principle models to predict an effective diffusion coefficient from the pore structure extracted from the images. The analysis also provides a rationale for why effective medium models do not predict accurate properties for this system. These alternate modeling approaches and their limitations as applied to our system are discussed below. The effective medium approach provides formulas that predict effective diffusivity coefficients from the porosity and tortuosity of the system. These formulas have a general form Deff = Dg(ε)/f(τ), where ε and τ are the porosity and tortuosity, respectively, and the functions “g” and “f ” are simple functions, typically low order polynomials. In its simple form Deff = Dε/τ. As a thought exercise, we calculated ε and τ from the image analysis results from Table 3. The resultant effective diffusion coefficients are 2 orders of magnitude larger than the experimental data. For example with the use of a simple porosity model and the data generated in Table 3 the estimated effective diffusivity coefficient is ∼0.36D for device S1_P2. The actual effective diffusivity coefficient that fits the release data is reported in Table 2 (∼0.004D). Similar disparate discrepancies are observed for all devices fabricated in this work shown in Table 1. These results confirm that such simple formulas cannot describe our system. Since said formulas are derived from first principle analysis it is important to rationalize why they fail in our case. The derivations of such formulas are based on a characteristic microscale that divides the domains into subdomains (i.e., into a porous medium).

Figure 20. A comparison of microCT and nanoCT for an implant sample (S6_P2). Left: MicroCT acquired at 40× magnification with a pixel size of 690 nm. microCT does not have the spatial resolution for separating particles at the submicrometer scale. After erosion−dilation to separate the API domain into clusters, each individual cluster will still contain multiple API particles. Right: NanoCT acquired with a pixel size of 64 nm. The same sample was analyzed with both techniques, albeit not at the exact same location in the device. NanoCT has the capability to resolve the actual particle sizes, while achieving a contrast difference between the drug particles, polymer, and nanocavities. The delineated particle clusters are shown as connected components in microCT agglomerates, even though they actually contain multiple API particles as revealed by the actual connectivity from nanoCT. 15342


Article


From a speed perspective microCT is of great benefit. It can acquire extensive sample data sets of as-fabricated implants in only a few hours. Another benefit is that microCT can effectively analyze both macro-scale features (e.g., whole devices, global defects) as well as microscale features (e.g., drug domain information), which are both of interest to pharmaceutical formulation teams. This great flexibility makes microCT useful not only for the study of extended-release devices, but also enables application to others systems of importance in the pharmaceutical industry, which in turn is a major driver toward adoption of the technique in industry. Comparatively, the widespread application of 3D nanoscale imaging toward pharmaceutical applications is still in its infancy. As these techniques continue maturing, more advance multiscale models can be formulated.

With the use of nanoCT it is demonstrated that our system is composed of an additional layer of domains in the submicrometer scale. These domains are the main resistance to diffusivity, and not treated in the effective medium model. First-principles modeling makes no modeling assumptions on the underlying pore structure, and the effective diffusivity coefficient is determined using random walk simulations of the segmented 3D image and perform detailed simulations of diffusion.18,24,28 First the porous topology is extracted from 3D XRCT images as described in section 4.2.1. The random walk simulations are simulated in the random structure and the spread of the particle population at large times is associated with the effective diffusivity coefficient (see refs 18, 24, and 28 for details of the implementation of this methodology). This approach has the advantage that relies on first principle calculations, and if the topology is extracted correctly, the extracted effective diffusivity should match the experimental one. We implemented this method, and these results are shown in Table 5. Again, when compared with the values of the

6. CONCLUSIONS We introduce a minimalistic modeling approach for which the effective release coefficient of a long-acting implant is estimated from drug domain sizes and coordination number distributions extracted from XRCT microstructural images. As for any simplified theory the method does not have universal validity and should be applied in the range where it was calibrated. This simple model provides a physical interpretation of parameters that can be extracted from an XRCT imaging, which can be used to assess drug release changes due to changes in implant microstructure (e.g., the role of drug microdomain size and connectivity on drug release). A detailed theoretical analysis justifies the model and supports some of the most salient model assumptions using image analysis. The model explains very well the experimental data consisting of implants fabricated by two different equipment setups, different drug loadings, and different process conditions (including varying temperatures and shear rates). The case studies demonstrate that this modeling strategy can provide practical guidance during process development, reducing some of the long lead time needed to generate release curves of extended-release formulations empirically. Furthermore, the model gives insight to the development engineers regarding the microstructural changes that have occurred, which is not available from an empirical release curve. Although the outlined strategy is very practical, to fully understand the pore structure, higher resolution is needed. We introduced for the first time nanoCT imaging of pharmaceutical implants. Using nanoCT we demonstrate a 10-fold improvement in resolution, revealing finer intrapore structures not captured by microCT. The new features indicate that although the minimalistic model based on microCT can provide good insight into the underlying factors that control the release and predict very well the release behavior, we should be careful in extrapolating the calibrated model to predict extreme conditions, as the submicrometer structure detected by nanoCT may change in a nontrivial way for very extreme process changes. Future work conducted within Merck and with external collaborator will consider more advanced transport modeling and image-analysis techniques, specifically to enable a direct application of pore-scale models toward release modeling by quantification of 3D nanoscale data.

Table 5. A Comparison of the Effective Diffusivity Coefficients Determined by a Random Walk Method from a Direct Application of the XRCT Microstructure, versus Experimental Release Data sample

Deff/D (random walk)

Deff/D (experimenta fit)

S1_P1 S2_P2 S1_P2 S2_P2

0.28 0.37 0.28 0.29

0.0024 0.0080 0.0004 0.0015

experimental fit the theoretical values are 2 orders of magnitude higher. Again, the discrepancies between theoretical and experimental values are due to the finer nanoscale structures that are needed to make these first principle calculations precise. This hinders the applicability of microscale models to our system. To produce a modeling methodology that can predict the effective diffusivity coefficient from XRCT images, we artificially reduced the interparticle connection diameter to the geodesics tracks to explore how the effective diffusivity coefficient increases with reduction of neck diameter. These numerical experiments resulted in pore diameters of about 30 nm that is in the same order of magnitude of the pore diameter of about 60 nm, providing further evidence that submicrometer resolution is needed to use this method effectively. This was also confirmed experimentally by the use of NanoCT. The improvements in image resolution afforded by NanoCT demonstrate that although the minimalistic model based on microCT can estimate the release behavior and gives good insight into which process parameters have a major impact on changes in drug release kinetics, we should be careful in extrapolating the calibrated model to predict extreme conditions, as the nanostructures detected by the nanoCT may change in a nontrivial way for very extreme process changes. For example, similar pore size distribution at the microCT resolution may have radical intrapore differences at the nanoCT resolution. This means that future model enhancements should include nanoCT with finer resolution to further widen the operating window of a calibrated model. With higher resolution images, more advanced image analysis and modeling techniques will be required, which are out of the scope of this article but are part of ongoing work. 15343


Article


■

APPENDIX A: MACROSCOPIC DIFFUSIONAL MODEL The geometry for the macroscopic model is shown in Figure 4. It is assumed that a diffusion process is symmetric in cylindrical coordinates. ∂C ∂ i ∂C yz 1 ∂(Deff C) zz + = jjjDeff ∂t ∂r k ∂r { r ∂r

(A1)

C(hmov ) = Cs

(A2)

C(R dev) = 0

(A3)

dWacc = JA dt

(A4)

ij D yz ∂C dh = jjj eff zzz j Cs z ∂r dt k {

of samples; Leon Farber for learnings and support on XRCT imaging; Zhen Liu for his work to expand imaging and modeling capabilities; Randolph Crawford for assisting in performing measurements of nanoCT images; Pavithra Sundararajan, Zhen Yang, Yash Kapoor, Binfeng Xia, Jingtao Zhang, Jochen Schoell, Richard J. Varsolona, and Eric L. Margelefsky for their interest in and discussion about the work; Timothy Rhodes for organizing the symposium that served as the starting point for these collaborations, as well as Wei Xu, Stephanie E. Barrett, Matthew Lamm, Ansu Bagchi, Belma Dogdas, Marian Gindy, and Anthony Leone for their leadership support to advance imaging and modeling capabilities in pharmaceutical formulation development. We would also like to thank Zeiss for enabling the nanoCT analysis, specifically Steve Kelly for data acquisition and Steve Neptun for coordinating evaluation activities.

■

(A5)

Numerical Solution

The moving boundary problem was mapped to a fixed domain ξ ∈ [0,1] using the following bilinear mapping:45 ξ=

r − hmov R dev − hmov

(A6)

This mapping fixes the interface at ξ = 0, and transforms the changing domain r ∈ [hmov, Rdev] into the transformed fixed domain ξ ∈ [0,1]. Equations A1−A5 are rewritten in terms of the transformed domain using the transformation relation eq 6. The time derivative in eq A1 is modified as follows to account for the effect of the moving nodes in the physical domain: ∂C(ξ , t ) dξ(t ) ∂C(ξ , t ) ∂C(r , t ) + = ∂ξ ∂t ∂t dt

(A7)

The transformed equation in the fixed domain are discretized using second order finite difference and integrated numerically using Runge−Kutta with variable time control. The initialization of the integration of the ODE system after discretization is very important for numerical stability. The following initialization is utilized. At time zero, there is no diffusion domain hmov = Rdev. The simulations start at a short time Δt in which the front is now, hmov < Rdev. With this assumption the initial profile is a linear profile. From this point on the discretized model is integrated to generated concentration profiles, front position, and % of drug released as a function of time.

■

REFERENCES

(1) Maniruzzaman, M.; Boateng, J. S.; Snowden, M. J.; Douroumis, D. A review of hot-melt extrusion: process technology to pharmaceutical products. ISRN Pharm. 2012, 2012, 436763. (2) Baghel, S.; Cathcart, H.; O’Reilly, N. J. Polymeric Amorphous Solid Dispersions: A Review of Amorphization, Crystallization, Stabilization, Solid-State Characterization, and Aqueous Solubilization of Biopharmaceutical Classification System Class II Drugs. J. Pharm. Sci. 2016, 105 (9), 2527−44. (3) Musiime, S.; Muhairwe, F.; Rutagengwa, A.; Mutimura, E.; Anastos, K.; Hoover, D. R.; Qiuhu, S.; Munyazesa, E.; Emile, I.; Uwineza, A.; Cowan, E. Adherence to Highly Active Antiretroviral Treatment in HIV-Infected Rwandan Women. PLoS One 2011, 6 (11), 1−6. (4) Ramirez Garcia, P.; Cote, J. K. Factors affecting adherence to antiretroviral therapy in people living with HIV/AIDS. Journal of the Association of Nurses in AIDS Care: JANAC 2003, 14 (4), 37−45. (5) Whetten, K.; Reif, S.; Whetten, R.; Murphy-McMillan, L. K. Trauma, mental health, distrust, and stigma among HIV-Positive persons: Implications for effective care. Psychosom. Med. 2008, 70 (5), 531−538. (6) Siegel, R. A.; Langert, R. A new Monte Carlo approach to diffusion in constricted porous geometries. J. Colloid Interface Sci. 1986, 109, 426−440. (7) Siegel, R. A.; Langer, R. Mechanistic studies of macromolecular drug release from macroporous polymer. II. Models for the slow kinetic of slow release. J. Controlled Release 1990, 14, 153−167. (8) Costa, P.; Lobo, M. S. Modeling and comparison of dissolution profiles. Eur. J. Pharm. Sci. 2001, 13, 123−133. (9) Dash, S.; Murthy, P. N.; Nath, L.; Chowdhury, P. Kinetic model on drug release from controlled drug delivery systems. Acta Poloniae Pharmaceutica-Drug Reseach 2010, 67, 217−223. (10) Siepmann, J.; Siepmann, F. Modeling of diffusion controlled drug delivery. J. Controlled Release 2012, 161, 351−362. (11) Peppas, N. A.; Narasimhan, B. Mathematical models in drug delivery: How modeling has shaped the way we design new drug delivery systems. J. Controlled Release 2014, 190, 75−81. (12) Ferrero, C.; Massuelle, D.; Jeannerat, D.; Doelker, E. Towards elucidation of the drug release mechanism from compressed hydrophilic matrices made of cellulose ethers. III. Critical use of thermodynamic parameters of activation for modeling the water penetration and drug release processes. J. Controlled Release 2013, 170, 175−182. (13) Herrmann, S.; Winter, G.; Mohl, S.; Siepmann, F.; Siepmann, J. Mechanisms controlling protein release from lipidic implants: Effects of PEG addition. J. Controlled Release 2007, 118, 161−168. (14) Higuchi, T. Rate of release of medicaments from ointment bases containing drugs in suspension. J. Pharm. Sci. 1961, 50, 874− 875.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Daniel Skomski: 0000-0002-8401-5626 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge the support from a number of co-workers at Merck research laboratories: Wei Xu for her tremendous advisorial and managerial support to pursue this project as well as many valuable teachings; Stephanie E. Barrett for her stewardship and encouragement; Hanmi Xi for her supervision; Carol Freddo for microtoming 15344


Article

Industrial & Engineering Chemistry Research (15) Siepmann, J.; Peppas, N. A. Higuchi equation: derivation, applications, use and misuse. Int. J. Pharm. 2011, 418, 6−12. (16) Paul, D. R. Elaborations on the Higuchi model for drug delivery. Int. J. Pharm. 2011, 418, 13−17. (17) Xiong, Q.; Baychev, T. G.; Jivkov, A. P. Review of pore network modelling of porous media: experimental characterisations, network constructions and applications to reactive transport. J. Contam. Hydrol. 2016, 192, 101−117. (18) Blunt, M. J.; Bijeljic, B.; Dong, H.; Gharbi, O.; Iglauer, S.; Mostaghimi, P.; Paluszny, A.; Pentland, C. Pore-scale imaging and modelling. Adv. Water Resour. 2013, 51, 197−216. (19) Schlüter, S.; Sheppard, A.; Brown, K.; Wildenschild, D. Image processing of multiphase images obtained via X-ray microtomography: a review. Water Resour. Res. 2014, 50 (4), 3615−3639. (20) Al-Kharusi, A. S.; Blunt, M. J. Network extraction from sandstone and carbonate pore space images. J. Pet. Sci. Eng. 2007, 56 (4), 219−231. (21) Tompson, A. F. B.; Gelhar, L. G. Numerical simulation of solute transport in three-dimensional, randomly heterogeneous porous media. Water Resour. Res. 1990, 26, 2541−2562. (22) Wu, R.; Kharaghani, A.; Tsotsas, E. Two phase flow with capillary valve effect in porous media. Chem. Eng. Sci. 2016, 139, 241− 248. (23) Blunt, M. J.; Bijeljic, B.; Dong, H.; Gharbi, O.; Iglauer, S.; Mostaghimi, P.; Paluszny, A.; Pentland, C. Pore-scale imaging and modeling. Adv. Water Resour. 2013, 51, 197−216. (24) Xiong, Q.; Jivkov, A. P.; Yates, J. R. Discrete modelling of contaminant diffusion in porous media with sorption. Microporous Mesoporous Mater. 2014, 185, 51−60. (25) Xiong, Q.; Baychev, T. G.; Jivkov, A. P. Review of pore network modeling of porous media: experimental characterizations, network constructions and application to reactive transport. J. Contam. Hydrol. 2016, 192, 101−117. (26) Bekri, S.; Thovert, J. F.; Adler, P. M. Dissolution of porous media. Chem. Eng. Sci. 1995, 50, 2765−2791. (27) Metzger, T.; Irawan, A.; Tsotsas, E. Influence of Pore Structure on Drying Kinetics: A Pore Network Study. AIChE J. 2007, 53, 3029−3041. (28) Liu, M.; Mostaghimi, P. High-resolution pore-scale simulation of dissolution in porous media. Chem. Eng. Sci. 2017, 161, 360−369. (29) Peng, S.; Hassan, A.; Loucks, R. G. Permeability estimation based on thin-section image analysis on 2D flow modeling in graindominated carbonates. Mar. Pet. Geol. 2016, 77, 763−775. (30) Okabe, H.; Blunt, M. J. Prediction of permeability for porous media reconstructed using multiple-point statistics. Phys. Rev. E 2004, 70 (6), No. 066135, DOI: 10.1103/PhysRevE.70.066135. (31) Piri, M.; Blunt, M. J. Three-dimensional mixed-wet random pore-scale network modeling of two-and three-phase flow in porous media. I. Model description. Phys. Rev. E 2005, 71 (2), No. 026301, DOI: 10.1103/PhysRevE.71.026301. (32) Jiang, Z.; Wu, K.; Couples, G.; Van Dijke, M. I. J.; Sorbie, K. S.; Ma, J. Efficient extraction of networks from three-dimensional porous media. Water Resour. Res. 2007, 43 (12), No. 5780, DOI: 10.1029/ 2006WR005780. (33) Dong, Hu; Blunt, M. J. Pore-network extraction from microcomputerized-tomography images. Phys. Rev. E 2009, 80 (3), No. 036307, DOI: 10.1103/PhysRevE.80.036307. (34) Bhattad, P.; Willson, C. S.; Thompson, K. E. Effect of network structure on characterization and flow modeling using X-ray microtomography images of granular and fibrous porous media. Transp. Porous Media 2011, 90 (2), 363−391. (35) Jivkov, A. P.; Hollis, C.; Etiese, F.; McDonald, S. A.; Withers, P. J. A novel architecture for pore network modelling with applications to permeability of porous media. J. Hydrol. 2013, 486, 246−258. (36) Tahmasebi, P.; Hezarkhani, A.; Sahimi, M. Multiple-point geostatistical modeling based on the cross-correlation functions. Comput. Geosci. 2012, 16, 779−797.

(37) Babaei, M.; Joekar-Niasar, V. A transport phase diagram for pore-level correlated porous media. Adv. Water Resour. 2016, 92, 23− 29. (38) Liu, M.; Mostaghimi, P. Characterisation of reactive transport in pore-scale correlated porous media. Chem. Eng. Sci. 2017, 173, 121−130. (39) Auriault, J. L. Effective macrospocpic description for heat conduction in periodic composites. Int. J. Heat Mass Transfer 1983, 26, 861−869. (40) Hassani, B.; Hinton, E. A review of homogenization and topology optimization, I-homogenization theory for media with periodic structure. Comput. Struct. 1998, 69, 707−717. (41) Otsu, N. A threshold selection method from gray-level histograms. Automatica 1975, 11 (285−296), 23−27. (42) Haralick, R. M.; Sternberg, S. R.; Zhuang, X. Image analysis using mathematical morphology. IEEE transactions on pattern analysis and machine intelligence 1987, PAMI-9 (4), 532−550. (43) Tenenbaum, J. B.; De Silva, V.; Langford, J. C. A global geometric framework for nonlinear dimensionality reduction. Science 2000, 290 (5500), 2319−2323. (44) Irizarry, R. LARES: An Artificial Chemical Process Approach for Optimization. Evolutionary Computation Journal 2004, 12 (4), 435−460. (45) Derby, J. J.; Brown, R. A. A fully implicit method for simulation of the one-dimensional solidification of a binary alloy. Chem. Eng. Sci. 1986, 41, 37−36.

15345


Theoretical Modeling and Mechanism of Drug Release from Long

Recommend Documents