Impact of Copolymer Architecture on Demicellization and Cargo Release via Head-to-Tail Depolymerization of Hydrophobic Moieties

1. Introduction

Hydrophobic drugs and nucleic acid vectors used in therapeutics face significant challenges, including low stability, in vivo degradation, and difficulties in reaching the target site [1]. Consequently, a robust delivery system is essential to transport these agents safely and efficiently. In recent years, polymeric micelles have emerged as highly effective cargo carriers, owing to their biocompatibility, facile preparation, tunable size and shape, prolonged stability in blood circulation, and enhanced in vivo retention at the target site [2,3,4,5,6]. Typically, polymeric micelles consist of pre-synthesized diblock copolymer chains, comprising a hydrophilic and a hydrophobic block, or a mixture of double hydrophilic diblock copolymers with neutral and oppositely charged blocks, that become hydrophobic upon complexation [7,8]. The hydrophobic segments form the micelle core, while the hydrophilic ones create the surrounding corona. These micelles are derived from dilute polymer concentrations exceeding the critical micelle concentration (CMC). Recently, innovative assembly models (PISA and PIESA) have allowed micellization at high monomer concentrations, combining diblock copolymer synthesis and micellization in a single process [9,10]. The coassembly of cargo molecules with copolymer chains results in the encapsulation of hydrophobic cargo within the micelle core due to the hydrophobic interactions with the solvent.

Smart polymeric micelles, responsive to the microenvironment of the target site, can be designed by incorporating stimuli-responsive switches either assembled into micelles or covalently grafted onto copolymer chains. These switches, responsive to pH, light, heat, redox reactions, and enzymes, facilitate efficient drug delivery [11,12,13]. However, challenges may arise in achieving complete cargo release due to the incomplete dissociation of the micelle. Polymeric micelle carriers with degradable copolymer chains have gained attention for offering controlled yet complete release of cargo molecules and preventing the accumulation of micelles at the target point [14,15,16].

In a noteworthy study Kim et. al. designed amphiphiles based on poly(benzyl-ether) capable of selective demicellization through head-to-tail depolymerization triggered by hydrophobic fluolid molecules [17]. A labile end-capping unit attached to the end of hydrophobic block is buried inside the core of the micelle, preventing detachment of the labile units and rendering the demicellization signal specific. Fluolid molecules trigger the detachment of the end-capping unit in the hydrophobic block, after which the entire chain spontaneously and continuously depolymerizes in a head-to-tail manner without the need for additional stimuli. Loading doxorubicin inside the micelles led to their molecular-level degradation, resulting in controlled and complete release of the cargo molecules.

The depolymerization-induced disassembly mechanism has primarily been explored for amphiphiles with a linear chain architecture carrying one labile end cap unit. The design of amphiphiles with a different architecture, bearing more labile end cap units, can further enhance the control of depolymerization and the cargo release, which is useful for biological and environmental applications. Due to the lack of such studies, we conducted coarse-grained molecular dynamics simulations of micelles formed by a) Linear AB diblock copolymer with one hydrophilic (A type) and one hydrophobic (B type) block, featuring one end cap unit. b) Miktoarm star copolymer A(B)₂ with one hydrophilic and two hydrophobic branches, equipped with two end cap units. c) Miktoarm star copolymer A(B)₃ with three hydrophobic branches and three end cap units. We calculated the mass distributions of micelles formed by different copolymers before and after cargo loading, the kinetics of the head-to-tail depolymerization of hydrophobic moieties and the kinetics of cargo release triggered by hydrophobic small molecules.

2. Model

2.1. Coarse-Grained Molecular Dynamics Simulation Details

Amphiphilic AB linear diblock copolymers, miktoarm star copolymers, and cargo molecules are modeled as bead-spring chains comprising Lennard-Jones beads with a diameter σ. In all simulations, the copolymer chains collectively consist of 30 A-type hydrophilic and 30 B-type hydrophobic beads. In linear AB diblock copolymers, denoted as A₃₀B₃₀, the hydrophilic and hydrophobic beads are distributed in one block each. In miktoarm star copolymers with 3 and 4 branches, denoted as A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃, the 30 hydrophobic beads are distributed in 2 or 3 branches, respectively. The free end bead of the hydrophobic block or branch is designated as an end cap unit. The cargo molecules are linear chains containing 3 C-type hydrophobic beads, denoted as C₃. The trigger molecule consists of a single T-type hydrophobic bead with a diameter σ. All copolymer chains, including those with linear and star architecture, as well as cargo and trigger molecules are illustrated in Figure 1.

In the current study, simulations were contacted for mixtures comprising a) A₃₀B₃₀ with C₃ b) A₃₀(B₁₅)₂ with C₃ and c) A₃₀(B₁₀)₃ with C₃, collectively containing 1000 copolymer chains and 2000 or 4000 cargo molecules. In all cases, the total concentration of beads (including copolymer and cargo beads) was maintained at [Φ]=0.12 where most micelles are formed.

Excluded volume, bead-bead interactions were considered to mimic the macroscopic solvent conditions. These interactions were calculated using a truncated Lennard-Jones potential [8,18]:

$$U_{LJ}\left(r_{ij}\right)=\left\{\begin{array}{c}4\epsilon \left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6-{\left(\frac{\sigma }{r_{cij}}\right)}^{12}+{\left(\frac{\sigma }{r_{cij}}\right)}^6\right],r_{ij}\le r_{cij}\\ 0,r_{ij}>r_{cij}\end{array}\right.$$

(1)

where ε is the well-depth and r_cij is the cutoff radius. Different beads were connected with finitely extensible nonlinear elastic bonds (FENE). The FENE potential is expressed as [8]

$$U_{Bond}\left(r_{ij}\right)=\left\{\begin{array}{c}-0.5k{R}_0^{2}ln\left[1-{\left(\frac{r_{ij}}{R_0}\right)}^2\right],r_{ij}\le R_0\\ \infty ,r_{ij}>R_0\end{array}\right.$$

(2)

where r_ij is the distance between beads i and j $,k=25\epsilon /{\mathsf{\sigma }}^2$ and R₀ is the maximum extension of the bond [19] ($R_0=1.5\mathsf{\sigma })$. The solvent molecules are implicitly considered.

Molecular dynamics simulations with a Langevin thermostat were conducted in a cubic box with periodic boundary conditions, using the open-source massive parallel simulator LAMMPS [20]. The reduced temperature of the simulation T* was set to T*=k_BT/ε=1.8 (ε=1) corresponding to bad solvents conditions [8]. Different cutoff distances and epsilon parameters in the Lennard-Jones potential were employed to describe the interactions between copolymer units. The B-B, B-C, and C-C interactions were considered attractive. Specifically, the B-B interactions consistently corresponded to T*=1.8 while the B-C and C-C interactions were even more attractive, corresponding to lower temperatures (T*=1.8, 1.4 and T*=1.8, 1.6, 1.5 and 1.4 respectively) with a cutoff distance r_cij of 2.5σ (ε_BB=1, ε_ij= T*/1.8, for i, j≠B). Conversely, the A-A, A-B, and A-C interactions were considered repulsive with cutoff distance r_cij of 2^1/6σ. In the latter case, the Lennard Jones potential is shifted by ε. For simplicity, all types of beads were assumed to have the same mass (m= m_i =1) and diameter (σ=1).

Following the Stillinger criterion [21] two moieties were considered to reside to the same micelle if the distance between any two nonbonded hydrophobic beads, belonging to different chains, were found to be within 1.5σ.

We performed 2 million timesteps with an integration step $t=0.006\mathsf{\tau }$ ($\mathsf{\tau }=\sqrt{\frac{{\mathrm{m}\sigma }^2}{\epsilon }})$, setting all cutoff radii equal to r_cij=2^1/6σ to eliminate any bias introduced from the initial conformation. Subsequently, the system was allowed to equilibrate for 30 million timesteps. The simulation was then extended to 300 to 600 million timesteps. The duration of the simulation was evaluated by calculating the tracer autocorrelation function [7]:

$$C\left(t\right)=\frac{〈N\left(t_0+t\right)N\left(t_0\right)〉-{〈N\left(t_0\right)〉}^2}{〈N^2\left(t_0\right)〉-{〈N\left(t_0\right)〉}^2},$$

(3)

where N(t) is the number of molecules in the micelle in which the copolymer resides at time t. We considered all copolymers as tracers, with each time step serving as the time origin t₀. The characteristic relaxation time t_relax is defined as the time required for C(t) to reach the value [7,8] of 1/e = 0.37. Each simulation was conducted for a minimum of 10t_relax to ensure 10 independent conformations. The properties of interest were calculated as averages from 4000 snapshots using the block average method with ten blocks.

2.2. Depolymerization and Cluster Analysis

The stochastic depolymerization of the hydrophobic monomers of amphiphilic copolymers as depicted in Figure 2, utilized the “bond /react” functionality of LAMMPS. This feature allows the breaking of bonds between beads based on distance-dependent probabilistic criteria. Topology changes are specified in pre- and post-reaction molecule templates. Two distinct depolymerization steps were undertaken: (a) the breaking of bonds between the hydrophobic end cap and B-type beads, triggered by the presence of T-type molecule within R_cutoff = 1σ from the end cap bead, with a predetermined reaction probability (RP) (Figure 2a). (b) The spontaneous propagation of depolymerization of hydrophobic beads with a predetermined RP occurring without any external stimuli (Figure 2b).

We utilized graph clustering analysis to analyze the simulation data. Initially, we identified micelles employing the data clustering algorithm DBSCAN implemented in the Python library Sklearn [22], with a maximum allowable neighborhood radius of 1.5σ. For a point (bead) to be classified as a core point, a minimum of two points (including the point itself) must be within the neighborhood. A precomputed neighbor sparse array was utilized as input for the DBSCAN algorithm. To construct this array, the KDTree neighbor data structure from the Python library SciPy [23] was employed, specifically the Sparse Distance Matrix algorithm with a maximum distance of 1.5σ between two points. It’s worth noting that the distance matrix algorithm disregards points with a distance greater than the maximum distance parameter. The KDTree neighbor data structure takes periodic boundary conditions into account, ensuring that the clustering analysis includes the periodic images.

For the identification of polymer chains, we employed the Python NetworkX library [24]. In this context, beads were represented as nodes, and bonds were represented as edges. From the graph created by this library, we extracted the polymer chains using the “Connected/Components” algorithm. This algorithm generates connected components from a graph, corresponding to bead spring chains in our case. Subsequently, the polymer chains were assigned to micelles based on the previous steps. To compute properties such as the radius of gyration of the micelles (core, corona, and total) and the shape anisotropy parameter κ², we utilized the outbox coordinates. To achieve this, micelles that split due to the periodic conditions (inbox coordinates) were identified and consolidated using the data clustering algorithm DBSCAN, KDTree neighbor data structure, and Sparse Distance Matrix algorithm without periodic conditions.

3. Results and Discussion

3.1. Micellization and Cargo Encapsulation

The micellization of amphiphilic linear A₃₀B₃₀, miktoarm star A₃₀(B₁₅)₂, and A₃₀(B₁₀)₃ copolymers was studied prior to cargo encapsulation at T*=1.8 and [Φ]=0.12. We calculated the aggregate mass distribution, the gyration radius, and the shape anisotropy of micelles. In Figure 3, the mass distribution of the aggregates is depicted. Peak values corresponding to the preferential aggregation number N_p are higher for linear copolymers and decrease as the number of arms in miktoarm copolymers increases, aligning with previous simulation results [25]. This behavior can be attributed to varying interfacial energy per chain required from these copolymers during the aggregation process [26]. In the case of linear copolymers, the long hydrophobic arm may extend far from the connection point with the hydrophilic block, reducing the interface with the solvent and, consequently, the enthalpic interactions [27]. Additionally, the extension of the large hydrophobic block increases the interface exposed to other copolymers’ hydrophobic arms, enhancing the probability of interaction reducing the free energy. For miktoarm copolymers with two or three shorter arms, many beads are confined near the connection point with the hydrophilic arm. This leads to a larger interface with the solvent molecules, resulting in greater interfacial energy and a smaller aggregate. Furthermore, the concentration of beads near the star center hinders their interaction with other hydrophobic units, further reducing free energy.

The encapsulation of cargo molecules was examined for different interaction parameters between hydrophobic C-C and B-C beads. Simulations of mixtures containing 1000 copolymer chains with 2000 or 4000 cargo molecules C₃ were contacted at [Φ]=0.12, and the results are presented in Table S1. For linear A₃₀B₃₀ copolymers, an increase in the strength of attractions between C-C beads from T^*_C-C =1.8 to 1.4 had minimal effect on the percentage of encapsulated cargo molecules (41% to 43%). In contrast, a similar increase in the strength of attractions between B-C beads (T^*_B-B=1.8 and T^* _B-C = T^*_C-C=1.4) dramatically increased the percentage of encapsulated molecules to 83%, aligning with previous findings [28]. The copolymer architecture did not influence the encapsulation of cargo molecules, as the driving force for this process is the hydrophobicity between B and C beads. In a mixture of A₃₀B₃₀ with 4000 C₃ molecules and attractive interactions corresponding to T^*_B-B=1.8 and T^*_B-C = T^*_C-C=1.4, the percentage of encapsulated cargo molecules remained constant at 82%. Free cargo molecules were in dynamic equilibrium with the encapsulated ones. The loading capacity [17] (LC) of the micelle, defined as the w/w percentage of the cargo accommodated in the loaded micelle carrier for the two different cargo concentrations, was found to be 7.6% and 14.1%, respectively. Encapsulation of hydrophobic molecules in the micelle core altered the hydrophilic/hydrophobic balance, triggering the association of free copolymer chains into the micelles and resulting in a higher aggregation number for all copolymer architectures (Figure 3 and Table S2). The N_p of loaded micelles with LC=7.6% increased by approximately 16%, 30%, and 35% for the linear A₃₀B₃₀, miktoarm star A₃₀(B₁₅)₂, and A₃₀(B₁₀)₃ copolymer mixtures, respectively. This higher increase in N_p in miktoarm copolymers reflects the more efficient packing of B and C beads in the micelle core, facilitating the decrease in free energy. The autocorrelation function indicates that the micelles of A₃₀B₃₀ copolymers with LC=14.1% were kinetically frozen even after 600 million timesteps and were not further studied.

The mean squared radius of gyration <R_g²>, serves as a measure expressing the size of both ‘empty’ and loaded micelles. The calculation of <R_g²> focused solely on copolymer beads to maintain consistent molecular weight in micelles sharing the same aggregation number N. In the case of preferential micelles, <R_g²> was observed to be higher in loaded micelles by approximately 9%, 16%, and 16% for A₃₀B₃₀, A₃₀(B₁₅)₂, and A₃₀(B₁₀)₃ copolymers, respectively. This elevation can be attributed to the increased N_p of loaded micelles, with the disparity in N_p becoming more pronounced in miktoarm star copolymers. Figure 4 illustrates the results for <R_g²> for both types of micelles with the same N_p. For A₃₀B₃₀ copolymers, the <R_g²> is marginally higher in loaded micelles, a consequence of cargo molecules enlarging the micelle core and leading to a slight increase in <R_g²>. A similar trend is observed in micelles formed by A₃₀(B₁₅)₂, and A₃₀(B₁₀)₃ copolymers up to N=40 and 27, respectively. However, for higher N values the <R_g²> of ‘empty’ micelles surpasses that of loaded ones. The anisotropy shape parameter [27,29] κ² values, depicted in Figure 5 reveal that ‘empty’ micelles exhibit an elongated shape in contrast to the more spherical shape loaded micelles. This distinction arises from the improved packing of hydrophobic content in the micelle core, leading to higher <R_g²> values for ‘empty’ micelles compared to the loaded ones. Notably, A₃₀B₃₀ copolymers form very spherical micelles in both ‘empty’ and loaded states, while the differences in κ² values for A₃₀(B₁₅)₂, and A₃₀(B₁₀)₃ copolymers increase with the number of hydrophobic arms.

3.2. Kinetics of Depolymerization

The depolymerization and the subsequent degradation of loaded micelles were carried out using two distinct concentrations of hydrophobic trigger molecules. The first concentration aligns with the stoichiometry of the total end cap beads in A₃₀B₃₀, A₃₀(B₁₅)₂, and A₃₀(B₁₀)₃ copolymer solutions, accounting for 1000, 2000 and 3000 T-type molecules, respectively. The second concentration is maintained at a constant level, with 4000 T-type molecules applied to all copolymer solutions. These trigger molecules were randomly inserted into the simulation box, representing the final snapshot of the trajectory of copolymer comicellization involving cargo C₃ molecules. This configuration served as the initial state for the depolymerization process. The introduction of trigger molecules led to an adjustment in overall solution concentration, reaching up to [Φ]=12.6. The attractive interactions between T-T, B-T, and T-C beads consistently corresponded to T*=1.8 while the T-A interactions were considered repulsive.

When T-type molecules infiltrate the hydrophobic core of the micelle, they can initiate the rapture of bonds between the hydrophobic end cap and B-type beads if the distance of T-type molecule is within R_cutoff = 1σ from the end cap bead, governed by a predetermined reaction probability, RP_T. Subsequently, spontaneous head-to-tail depolymerization of hydrophobic B-type beads may occur, guided by a predetermined reaction probability, RP_B, and this process takes place without any external stimuli. Simulations were conducted to explore scenarios where the end cap bond-breaking occurs more easily, equivalently, or more difficultly than the B-B bond (with RP_B values of 10^-3, 10^-4 and RP_T values of 10^-2, 10^-3 and 10^-4).

The results depicting the fraction of end cap bonds broken with RP_T=10^-4 and RP_B=10^-3, are showcased in Figure 6 for both stoichiometric and constant concentrations of trigger molecules. Notably, the depolymerization of end cap beads in A₃₀(B₁₀)₃ copolymers at a stoichiometric concentration of trigger molecules surpasses the rates observed in A₃₀(B₁₅)₂ and A₃₀B₃₀ copolymers. This heightened efficacy is attributed to the doubled and tripled presence of T-type molecules in A₃₀(B₁₀)₃ mixtures compared to their counterparts in A₃₀(B₁₅)₂ and A₃₀B₃₀ mixtures, thereby increasing the probability of infiltration into the hydrophobic core. Furthermore, the elongated sphere scheme (depicted in Figure 5) of A₃₀(B₁₀)₃ micelles contribute to a larger core surface exposed to the solvent compared to other copolymers. This structural characteristic facilitates the penetration of trigger molecules, contributing to the observed faster depolymerization in A₃₀(B₁₀)₃ copolymers. Conversely, contrasting trends are noted in the depolymerization process under constant T-type molecules concentration (Figure 6b). The excess of trigger molecules is fourfold, twofold and 1.5 times higher than the stoichiometric case for A₃₀B₃₀, A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃ mixtures, respectively, leading to a quicker penetration of trigger molecules and depolymerization of linear copolymers followed by miktoarm star copolymers. Various factors, including micelle corrosion, also influence the penetration of trigger molecules into the micelle core and subsequently the breaking of the end cap bonds with B-type beads. Higher values of RP_T simply reduce the time for the completion of end cap removal without affecting the order of copolymer depolymerization in both stoichiometric and constant concentration cases (Figure S1). Similarly, an increase by x times in the T-type molecule concentrations in all copolymer mixtures results in faster end cap removal without altering the order of copolymer depolymerization (Figure S2).

The results of the depolymerization of hydrophobic beads, including both end cap and B-type beads, are presented in Figure 7 under constant and stoichiometric concentrations of trigger molecules. When the breaking of the end cap bond with the B-type bead is more challenging than the bond between B-type beads (RP_T=10^-4, RP_B=10^-3), the rate of depolymerization for hydrophobic beads is governed by the slower step. Consequently, the depolymerization trends among copolymer architectures align with those observed in the end cap bond break. Specifically, A₃₀B₃₀ exhibits faster depolymerization, followed by A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃ mixtures (Figure 7a) for the same reasons outlined earlier. However, at RP_T=10^-2 and RP_B=10^-3, where the breakage of the bond between B-type beads becomes the slower reaction step, the branching architecture significantly influences the depolymerization process. The head-to-tail depolymerization of B-type beads occurs simultaneously in all branches resulting in faster depolymerization in A₃₀(B₁₀)₃ with three branches, followed by A₃₀(B₁₅)₂ and A₃₀B₃₀ with two and one branches respectively. Notably, when 65% of the total bonds of hydrophobic beads are broken, the depolymerization rate of A₃₀(B₁₅)₂ copolymers surpasses that of A₃₀(B₁₀)₃ copolymers (Figure 7b). This is attributed to the higher excess of trigger molecules in A₃₀(B₁₅)₂ mixtures, facilitating the faster removal of end cap beads from the hydrophobic branch. Similarly, the depolymerization rate of A₃₀B₃₀ copolymers exceeds that of A₃₀(B₁₀)₃ and A₃₀(B₁₅)₂ copolymers when the fraction of total bonds breakage reaches 0.85 and 0.95 respectively. Beyond these values, the order of depolymerization rates reverts to the pattern observed at RP_T=10^-4 and RP_B=10^-3. In actuality, the reversal in the depolymerization order occurs in the case of RP_T=10^-4 and RP_B=10^-3, but the crossovers between the depolymerization rates are confined to a narrow range, specifically between 0.03 and 0.06. When RP_T=RP_B=10^-3, the depolymerization rate of hydrophobic beads for different architectures is determined by two opposite factors: the excess of trigger molecules, which increases as the number of hydrophobic branches or blocks decreases, and the simultaneous head-to-tail depolymerization of beads, which intensifies with an increasing number of branches. As depicted in Figure 7c, the crossover of depolymerization curves occurs within the range of 0.20 to 0.55. Beyond this fraction, A₃₀B₃₀ exhibits faster depolymerization, followed by A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃ mixtures.

In the case of depolymerization of hydrophobic beads under stoichiometric concentration of T-type molecules with the end cap beads (Figure 8), the rate of depolymerization is consistently faster in A₃₀(B₁₀)₃ followed by A₃₀(B₁₅)₂ and A₃₀B₃₀ mixtures. This can be attributed to the same factors explained earlier in the context of end cap depolymerization, coupled with the higher number of branches that are potentially subject to simultaneous head-to-tail depolymerization of B-type beads.

3.3. Demicellization

The depolymerization of hydrophobic moieties in a head-to-tail fashion results in copolymer chains that are either shorter or lack hydrophobic block or branches. The reduction in hydrophobic content within the copolymer corresponds to an increase in the critical micelle concentration. Consequently, these chains migrate into the solution, triggering the gradual demicellization of the mixtures. Figure 9, Figure 10 and Figure 11 and S3–S8 present the results regarding the micelle mass distribution at various fractions of bond breakage between hydrophobic beads (approximately 0, 0.15, 0.30, 0.45, 0.60, and 0.75) for stochiometric and constant concentrations of trigger molecules. Micelle mass distributions were computed from a single snapshot, and their normalization is based on the total micelle content at the specific snapshot. Copolymer chains devoid of B-type beads lose the ability to participate in micelles and are excluded from the count as single chains in the micelle distribution. Across all figures, it is evident that the aggregation number of micelles (N) progressively decreases with simulation time. This decline signifies the gradual loss of copolymer chains from micelles without micelle breaking.

In the case of A₃₀B₃₀ copolymers with stoichiometric mixtures of trigger molecules and the end cap beads (RP_T=10^-4 and RP_B=10^-3), Figure 9 illustrates that all chains leaving the micelles consist solely of A-type beads. This is attributed to the faster depolymerization of B-type beads in the B-type block of a copolymer chain due to the higher RP_B compared to RP_T. Conversely, in the case of single miktoarm copolymers, the chains exiting the micelles are composed of B-type moieties, characterized by a reduced content of B-type beads (Figure 10 and Figure 11). At a bond break fraction of 0.68, the aggregation number of A₃₀B₃₀ copolymer micelles remains substantial. However, in A₃₀(B₁₅)₂ and even more so in A₃₀(B₁₀)₃ copolymers, the aggregation numbers diminish at similar bond break fraction.

The same trends are observed in the micelles mass distribution for mixtures with a constant concentration of trigger molecules (Figures S3–S5). However, at the same bond break fraction, the aggregation numbers are noticeable smaller. When reaction probabilities for both reactions are equal (RP_T=RP_B=10^-3), micelles and single corroded copolymer chain are in equilibrium, but the aggregation number of micelles is smaller compared to the case with RP_T=10^-4 and RP_B=10^-3 at the same constant concentration. Increasing the reaction probability of end cap bonds (RP_T=10^-2 and RP_B=10^-3) leads to even smaller aggregation numbers for micelles of all architectures. Nonetheless, A₃₀B₃₀ micelles consistently exhibit higher aggregation numbers than the A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃ for the same bond break fraction (Figures S6–S8).

3.4. Kinetics of Cargo Release

In this section, the mobility of cargo molecules was investigated through the calculation of their mean squared displacement (MSD) as a function of time. The study utilized 90 cargo molecules residing within a micelle formed by A₃₀B₃₀ copolymers with N=50 for mixtures with RP_T=10^-4 and RP_B=10^-3, maintaining a constant concentration of trigger molecules. Results for MSD were averaged across the cargo tracers and are depicted in Figure 12. The exponent α corresponding to the scaling behavior of time was determined to be unity, suggesting normal diffusion in the release of cargo molecules. Notably, even before demicellization occurred in the tracer micelle, a substantial exchange of cargo molecules with different micelles and the solution was observed. Approximately, 18% of tracer molecules left the tracer micelle, while others were inserted, maintaining the number of encapsulation cargo molecules nearly constant. As the micelles underwent depolymerization and corrosion, the number of encapsulated cargo molecules decreased, while the number of molecules released into solution increased. The fraction ratio of released cargo molecules to the initially encapsulated ones in all micelles in the mixture is illustrated as a function of time in Figure 13 and Figure 9S. The observed noise in the curves is a result of released cargo molecules or others from the solution being inserted into other micelles. In general, the cargo release rates mirror the depolymerization rates of copolymers. Specifically, at a constant concentration of trigger molecules and RP_T=10^-4, RP_B=10^-3, the release of cargo molecules is faster in linear A₃₀B₃₀ followed by A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃ miktoarm copolymers. Conversely, for stoichiometric concentration of trigger molecules with the end cap beads, the release of cargo molecules is faster in A₃₀(B₁₀)₃ followed by A₃₀(B₁₅)₂ and linear A₃₀B₃₀ copolymers. This trend aligns with the decreasing hydrophobic content of B-type beads in micelles due to depolymerization and demicellization allowing fewer hydrophobic cargo molecules to be accommodated. However, the fraction of cargo molecules released to the solution is not directly proportional to the fraction of depolymerized B-type beads. For a 0.60 degradation of hydrophobic beads, the fraction of cargo release is approximately 0.34, 0.39 and 0.43 for A₃₀B₃₀, A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃, respectively at both stoichiometric and constant concentrations of trigger molecules with RP_T=10^-4 and RP_B=10^-3. This is attributed to the existence of micelles with higher N in A₃₀B₃₀ than in the miktoarm copolymer mixtures, allowing them to accommodate a larger number of cargo molecules (Figure 9, Figure 10 and Figure 11).

To quantitatively analyze the release of cargo molecules, the non-linear Korsmeyer-Peppas equation was employed [30,31]:

$$\frac{{{\rm M}}_t}{{{\rm M}}_{\mathrm{\infty }}}=k{t}^n$$

(4)

In this equation, M_t/Μ_∞ signifies the fraction of released cargo molecules, t denotes the release time, k is the kinetic constant (with dimensions of τ^-1), and n is the diffusional exponent for cargo release (dimensionless). The kinetic constant k primarily conveys information about the drug formulation’s characteristics, particularly those related to micelles, whereas n is related to the drug release mechanism. For 0.45<n<1 the release is non-Fickian, and both diffusion and micelle relaxation-corrosion contribute to the release mechanism [31] The analysis focused on the initial 60% of the release curve. Logarithmic plots for constant and stoichiometric trigger molecule concentrations (RP_T=10^-4 and RP_B=10^-3) are presented in Figure 14 and Figure S10, respectively. Notably, for constant trigger molecule concentration, n values exhibit a decreasing trend from 0.98 to 0.97 to 0.80 as the number of hydrophobic blocks or branches increases, whereas the constant k values manifest the opposite trend. For stoichiometric mixtures, the variation of n is non-monotonic, with values ranging from 0.77 for A₃₀B₃₀, 0.85 for A₃₀(B₁₅)₂, to 0.70 for A₃₀(B₁₀)₃, respectively, and similarly for k. Consequently, our findings strongly suggest that the release of cargo from degradable micelles follows a non-Fickian pattern, where both diffusion and corrosion mechanisms contribute to the release process. This observation aligns with the consistent n values obtained in a previous study involving the cationic drug release from four arm starblock copolymers where the release behavior was dominated by chain relaxation-erosion induced by ion exchange that was dependent on pH [32].

4. Conclusions

To elucidate the impact of copolymer architecture on demicellization and cargo release through head-to-tail depolymerization triggered by specific stimuli, we conducted comprehensive molecular dynamics simulations. The depolymerization reaction involved two distinct steps: (a) the rapture of bonds between the hydrophobic end cap and B-type beads, initiated by the presence of T-type molecules with a predefined reaction probability (RP_T); (b) depolymerization in the absence of external stimuli, with reaction probability (RP_B). Our investigation focused on linear A₃₀B₃₀ and miktoarm star A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃ copolymers, featuring one, two and three hydrophobic arms (and end caps) respectively. We explored these copolymers under constant and stoichiometric concentrations of trigger molecules relative to end cap beads, examining a range of RP_T and RP_B values. Our findings revealed that the preferential aggregation number of micelles (N_p) was higher for linear copolymers and decreased as the number of arms in miktoarm copolymers increased. Interestingly, copolymer architecture did not impact the loading capacity of micelles (LC). However, for LC=7.6%, loaded micelles exhibited an increase of approximately 16%, 30%, and 35% in N_p, and a 9%, 16%, and 16% rise in the radius of gyration for the A₃₀B₃₀, A₃₀(B₁₅)₂, and A₃₀(B₁₀)₃ copolymers respectively. The rate of depolymerization of hydrophobic beads under stoichiometric T-type molecules concentration was consistently faster in A₃₀(B₁₀)₃ followed by A₃₀(B₁₅)₂ and A₃₀B₃₀ mixtures. For constant trigger molecule concentrations, the depolymerization rate depended on two opposite factors: the excess of trigger molecules, which increased as the number of hydrophobic branches or blocks decreased, and simultaneous head-to-tail depolymerization, which intensified with an increasing number of branches. Initially, at low break fractions, A₃₀B₃₀ exhibited faster depolymerization, followed by A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃ mixtures. For higher break fractions, the order is reversed, depending on RP_T and RP_B values. Generally, cargo release rates mirrored copolymer depolymerization rates. However, the fraction of cargo molecules released was not directly proportional to the fraction of depolymerized B-type beads. For a 0.60 degradation of hydrophobic beads, the fraction of cargo release was approximately 0.34, 0.39 and 0.43 for A₃₀B₃₀, A₃₀(B₁₅)₂ and A₃₀(B₁₀)₃, respectively, at both stoichiometric and constant trigger molecule concentrations with RP_T=10^-4 and RP_B=10^-3. Utilizing the non-linear Korsmeyer-Peppas equation, our results strongly suggest that the release of cargo molecules from degradable micelles follows a non-Fickian pattern, where both diffusion and corrosion mechanisms contribute to the release process. These insights may inform the design of copolymers chains with more triggering groups, proving valuable in the development of delivery vehicles for bio and environmental applications.