Interpretable active learning meta-modeling for the association dynamics of telechelic polymers on colloidal particles

Understanding the association dynamics of multicomponent and multiscale materials composed of interacting colloidal particles with polymers is crucial in various biological and chemical processes such as wound healing, involving platelet-polymer plugs that assemble at the site of damage [1], drug delivery systems [2], food products [3], cosmetics [4], and waterborne coatings [5]. Waterborne coatings were introduced as a sustainable alternative to conventional coatings, which often contain hazardous solvents that pose risks to health and the environment. However, if not appropriately formulated, they may fail to meet desired performance and criteria, particularly rheological behavior during application as well as their shelf life [6]. Formulating effective waterborne coatings demands fine-tuning a range of competing design requirements. The impacts of various design choices on rheological properties are intertwined and multiscale, making the selection of new formulations largely empirical.

The components of waterborne coating fluids that most impact their rheology are colloidal latex particles and associative polymers, all suspended in water. The associative polymers are typically telechelic; they consist of a long hydrophilic backbone, commonly polyethylene oxide (PEO), capped at each end by relatively small hydrophobic alkyl groups (stickers). The stickers primarily associate with hydrophobic surfaces such as latex colloidal particles. Their attraction strength is determined by the number of methylene groups [7]. These associative chains, a.k.a rheology modifiers, can take four possible configurations in latex-polymer solutions, namely, (a) a bridge between neighboring particles, (b) a loop on one single particle, (c) a dangling chain with one end in the solution, and (d) a free chain in the solution. In addition, within a paint mixture, their association with other components of the waterborne coatings, such as pigment particles, was determined to have a negligible effect on rheology [8].

The typical associative polymers used in waterborne latex coating formulations are hydrophobically ethoxylated urethanes (HEURs). In the solution by themselves, these telechelic polymers usually aggregate into flowerlike micelles with a core of hydrophobic ends and a corona of stretched hydrophilic blocks. At moderate concentrations, stickers entering and leaving the adjacent flowers enable network relaxation. The most straightforward theory for such solutions was developed by Green and Tobolsky [9] and termed the “temporary network model.” In this relatively simple theory, the network relaxation is determined by a single relaxation time, denoted as $τ$⁠. When there is no flow, this relaxation time is exponentially dependent on the free energy cost, $ΔG$⁠, for the escape of the sticker from the micelle, given by the equation $τ= τ 0exp( Δ G/ k B T)$⁠. Annable et al. [10] observed that the relaxation times determined rheologically for the HEUR solutions agreed with the above equation, assuming that $ΔG$ increases linearly with the length of the sticker. The increase was estimated as approximately 1  $k BT$ per methylene group added to the sticker.

Unlike pure HEUR solutions, latex particle-HEUR solutions possess a broad distribution of relaxation times with nonterminal behavior. Wang and Larson [11] developed Brownian dynamics (BD) simulations of small particles with explicit, but over-simplified, dumbbell chains representing polymers to understand the structure and dynamics of latex particle-polymer suspensions. They demonstrated that even a simplified system shows at least four relaxation times, including (1) local relaxation of the HEUR loops, $τ H E U R$⁠; (2) translation of the loops along the particle surfaces, $τ s u r f$⁠; (3) relaxation of the HEUR bridges, $τ b r i d g e$⁠; and (4) relaxation of clusters of bridged particles, $τ c l$⁠. In general, $τ H E U R$ and $τ s u r f$ are very small (<1  $mus$⁠), and, therefore, these modes do not affect the experimentally detectable contributions to rheological behavior. These simulations facilitated the investigation of the interaction between polymers and particles. However, they can only predict stress relaxation for small particles with rapid diffusion, involving up to several dozen chains per particle. Simulating more realistic systems with larger particles (⁠ $∼$100 nm), each covered by hundreds of HEUR chains, is prohibitively expensive. Furthermore, a small system size limits the investigation of slow dynamics, such as relaxation of clusters of bridged particles and of the system-spanning network.

Thus, to mitigate the computational challenges posed by the numerous degrees of freedom, Hajizadeh et al. [12] developed a novel hybrid population balance-Brownian dynamics (Pop-BD) model. Unlike explicit BD simulations that track each individual chain, Pop-BD replaces the chains with dynamic bonds/potentials that are continuously added or removed based on particle motion. This approach employs probabilistic rules to estimate the likelihood of bridge or loop formation and rupture, depending on the current state of the system. Reducing the systems’ degrees of freedom, the Pop-BD model can simulate systems with hundreds of large colloidal particles and hundreds of implicit chains per particle. The ability to capture relaxation times related to chain relaxation is constrained in Pop-BD because it does not account for the presence of the explicit chains. However, this innovative approach is particularly useful for exploring large-scale systems and their larger relaxation times of relevance to experimentally measurable rheological responses of realistic systems, while being informed by the essential physics of the polymer-particle interactions from high resolution models.

However, accurately informing the Pop-BD model requires knowledge of bridge populations and loop-to-bridge and bridge-to-loop transition rates [represented by mean loop lifetime (MLL) and mean bridge lifetime (MBL) in this study]. These transition rates are functions of various factors, including but not limited to the gap between two particles, the sticker-free energy or attraction strength and the gap-dependent spring potential for chain stretch. Tripathi et al. [13] proposed a constitutive model to capture the dynamic behavior of telechelic polymer networks encompassing micelles and polymeric chains. Quantifying the shear-dependent behavior of pure HEUR solutions, they approximated the loop and bridge transition rates between micelles. These rates were later used as inputs for the Pop-BD model developed by Hajizadeh et al. [12] It should be noted that the simulation of particles using the gap-dependent transition rates required for the Pop-BD simulations has now been integrated into HOOMD-blue [14], a versatile particle simulation toolkit [15]. In their 2018 study, Rezvantalab and Larson [16] investigated transitioning from loops to bridges in bead-spring chains within a dilute regime, neglecting polymer-excluded volume effects. Additionally, they assumed a constant bridge-to-loop rate irrespective of the extent of chain stretching. Their investigation underscored the necessity of achieving resolution at the Kuhn step level to obtain reliable findings. Subsequently, Travitz et al. [17] employed BD simulations to measure the faster bridge-to-loop transition time. This analysis accounted for the excluded volume effects at the level of a Kuhn step. Furthermore, they applied a self-consistent field theory (SCFT) developed by Zhang et al. [18] to determine the equilibrium ratio of the number of bridges and loops and used detailed balance to thereby infer the slower loop-to-bridge transition times. These simulations treated colloid particles as flat surfaces, which is valid only in the small interparticle gap-to-colloid radius ratios, where the Derjaguin approximation can be used to account for the nonuniform gap. This model was not verified based on explicit-chain BD simulations with smaller colloid particles.

The main challenge in developing a precise predictive model of waterborne coating fluids is the fact the chains’ association dynamics and separation distances of the colloids unfold across substantially different time and length scales. These scales span from microseconds and $∼$1 nm for the association dynamics to $∼$100 nm and seconds for re-arrangement of the colloidal particle network. Consequently, a single simulation proves inadequate for encapsulating the extensive array of time and length scales inherent in macroscopic many-colloid simulations. Nevertheless, simulation of these dynamics at both the smaller and larger scales is essential for accurate prediction of the rheological behavior of these materials. Thus, the outstanding challenge necessitates developing an integrated and refined predictive model for the association dynamics that can reliably relay results from fine-grain simulations to coarser-grained simulations, such as Pop-BD.

Addressing this challenge has long been a need, but recent advances in computer modeling across multiple scales, combined with the advent of advanced machine learning (ML) techniques, now offer promising solutions [19,20]. Specifically, active learning meta-modeling can significantly enhance the efficiency and accelerate the development of precise meta-models or surrogate models of expensive-to-evaluate simulations. These ML methodologies commonly integrate Gaussian process (GP) models [21] with adaptive sampling strategies, efficiently exploring the design space and building accurate surrogate models. The sampling strategies involve iteratively selecting the most informative simulations or experiments for training and targeting areas where the model’s predictions are uncertain or need further understanding [22,23]. Furthermore, there is a particular interest in understanding how ML-based meta-model predictions are influenced by system parameters from the design space. This requires interpretability techniques such as feature importance analysis [24], partial dependence plots, and SHAP (Shapley additive explanations) values [25]. Employing interpretability techniques can leverage the understanding of the output of the ML meta-model and uncover underlying patterns and interactions within the system parameters. These approaches provide crucial guidance for material design and discovery [26] with valuable insights into the fundamental mechanisms of the system.

Therefore, we initially aimed to develop new fine-grained simulations that maintain accuracy without relying on the simplifying approximations commonly employed in the literature [27]. Although computationally demanding, these predictive models can predict the association dynamics of polymeric chains across various system parameters, including chain length and population, the radius and spacing of latex particles, and the attraction strength between stickers and particles. In what follows, after developing these detailed simulations, we implemented novel active learning meta-models. The active learning process is able to minimize the number of samples necessary for building an accurate surrogate model over the design space, as it iteratively selects the most informative samples based on the knowledge gained from previous ones. This study prioritized interpretability throughout the model development process to promote the comprehension of these meta-models. The developed meta-models, consequently, facilitate the integration of the microscopic simulations into macroscale models by passing information from the detailed and fine-grained BD simulations to coarse-grained models such as Pop-BD. The gained insights also provide a flexible tool for exploring and optimizing the design parameters of waterborne coatings.

The paper is structured as follows. In Subsection II A, we investigate the dynamics of chains at the high resolution of a single Kuhn step (⁠ $∼$1 nm) by performing BD simulations of a minimal model, which allowed us to capture key responses, including loop and bridge lifetimes as well as bridge populations, essential for informing the lower-resolution simulations such as Pop-BD. In Subsection II B, we define the design space and develop an intelligent search algorithm to gradually learn the model responses across the design space. In Sec. III, we analyzed our findings to explain the results and highlight the significance of different parameters on the model responses using a machine learning interpretability tool. Finally, the conclusion and future direction of this research are discussed in Sec. IV.

The simulated model features both colloidal particles and HEUR polymeric chains. The particle size varies from around 10 nm in theoretical studies [11,27] to approximately 60 nm in typical latex suspensions. As previously noted, HEURs are telechelic polymers characterized by PEO backbones, end-capped on both sides with hydrophobic alkyl groups (stickers). The molecular weight of the PEO is 137 g/mol per each Kuhn segment with an effective length (⁠ $b k$⁠) of approximately 1.1 nm [28]. Experimental studies [5] typically employ telechelic polymers with a molecular weight of approximately 30 000 g/mol, corresponding to around 200 Kuhn steps. The strength of attraction (⁠ $ε s$⁠) between the stickers and colloidal particles is quantified by the number of carbon groups in the hydrophobic end, each contributing to approximately 0.75  $k BT$ [6]. The range of $ε s$ varies from 7 to 10  $k BT$ in this study. We restricted the hard-core radius of colloidal particles, $R p$⁠, to between 11 and 33 nm and the number of Kuhn steps (⁠ $N k$⁠) from 10 to 30. This parameter range resulted in end-to-end distance to colloid radius ratios of approximately between 0.1 and 0.5, which aligns with the ratios used in experiments (⁠ $∼$0.23). However, it should be noted that the ranges of $ε s$⁠, $R p$⁠, and $N k$ are restricted due to the current computational limitations. The number of chains per particle was assumed to range from 10 to 100, with the upper bound approaching typical experimental values [29]. The ranges of these variables are summarized in Table I.

To carry out the expensive fine-scale BD simulations resolved at the level of the Kuhn step needed to obtain the data for the Pop-BD model, we develop a “minimal model” in a periodic simulation box containing both polymeric chains and a fixed particle positioned at the center of the box. This model is the smallest unit that can represent a system able to allow simulation of bridge and loop transition dynamics needed for the Pop-BD simulations, with the periodic images of the particle serving as the neighboring particles to which bridges can form. The simulation box is cubic with a dimension L that we vary from run to run to set the gap $H$⁠, as $L=H+2× R t r u e$⁠, in which $H$ is then the effective interparticle gap, defined as the minimum distance between the surfaces of the particle and its periodic image, and $R t r u e$ is the true radius of the particles defined in Eqs. (3) and (4). The chains are initialized randomly in the space between particles where they can form the four possible configurations (loop, bridge, dangling, and free) depending on the stickers’ distance from the particle(s) surfaces. A sticker is considered “stuck” to a colloidal particle when it is within $R c u t , − 1$ of the particle’s center. $R c u t , − 1$ is defined as the distance at which the interaction energy is stronger than 1 $k BT$⁠, as depicted in Fig. 1. A chain is in the “loop” configuration if both of its stickers are stuck to the same particle, or in the “bridge” configuration if the stickers are stuck to different particles. Similarly, the chain is dangling or free, if one or none, of the stickers are stuck to a particle, respectively. The chains continuously transition between these four configurations; however, a similar study [32] has shown that chains spend the majority of their time in either the loop or bridge configuration for the range of $ε s$ considered here.

The critical results extracted from the simulations include the MBL (or bridge-to-loop transition time), the MLL (or loop-to-bridge transition time), and the MBF. The bridge lifetime measures how long chains stay in a bridge configuration before one end-sticker escapes from its particle’s attraction potential (exceeding $−$1  $k BT$⁠) and attaches to another particle, forming a loop. Conversely, the loop lifetime refers to the duration chains remain in a loop configuration, with both stickers attached to a single particle, until one end-sticker escapes the potential well (exceeding $−$1  $k BT$⁠) and forms a bridge. Any intermediate dangling or free chains that form during the transition from loop to bridge or vice versa are assumed to be part of the preceding bridge or loop configuration. Thus, there are only two states: bridge and loop, with transitional time periods between the two being lumped into the state existing prior to completion of the transition to the second state, following the same procedure used in the work of Travitz and Larson [17]. Before computing the average lifetimes, bridge and loop lifetime samples should be collected over a sufficiently long simulation time to ensure adequate statistical sampling. Further details regarding the sampling are discussed in Appendix A 1.

In this subsection, we detail the development of active learning meta-models of the MLL, mean bridge lifetime (MBL), and bridge fractions (BFs) that are extracted from the detailed BD simulations of colloidal particles interacting with “rheology modifying” telechelic polymer chains discussed in Sec. II A. As mentioned there and represented in Table I, the design space explored in this study is represented by five system parameters, namely, the hard-core radius of particles $R p/ b k∈[10−30]$⁠, number of Kuhn steps $N k∈[10−30]$⁠, normalized gap $H/ N k b k∈[0.1−0.6]$⁠, number of polymer chains per particle $N c/ N p∈[10−100]$⁠, and Lennard-Jones attraction strength between the stickers and the particles, where $ε s/ k BT∈[7−10]$⁠. While the ranges of particle sizes, chain lengths, and numbers of chains per particle to be explored in Table I are lower (by a factor of around 10) than the commercially most common, these values chosen to allow faster simulations and serve to provide a proof of principle. The ranges can be expanded upon in future work.

The active learning modeling aims at learning the model responses ${ f 1, f 2, f 3}$ = ${$MLL, MBL, BF $}$ across the entire design space, while reducing the number of simulations or experiments needed. If we were to conduct conventional grid-search (a form of exhaustive and unguided search) to find responses across the design space, it would have been necessary to perform a minimum of 1080 simulations for considering three options for each $R p$ and $N k$⁠, six options for $H/ N k b k$⁠, five options for $N c/ N p$⁠, and four options for $ε s$⁠. Our active learning meta-modeling aims to significantly reduce this exhaustive set of simulations by using an intelligent search methodology. The overall process is illustrated in Fig. 2, with each section explained in more detail in what follows. After defining the design space, the process begins by fitting initial meta-models $f ^ 1, f ^ 2, f ^ 3$ for initial random points $X= { x 1 , … , x n } T$ within the design space, where each $x i$ is a set of values of the system parameters or features (i.e., five parameters defined above). The essence of the process lies in making informed decisions about the next samples to simulate or experiment providing the most informative knowledge. This is guided by an acquisition function tailored to the specific problem. Once the new samples or a batch of new samples are selected, the BD simulations (or experiments) are carried out to update the fitted meta-models. This iterative process continues until specific stopping criteria are satisfied.

A maximin distance method [33] was used to generate a diverse and well-distributed set of initial random points, reducing the risk of clustering within the design space. The algorithm generates 1000 points and selects an initial set with the maximum pairwise distance between points, where “distance” is measured numerically in terms of the differences in values of the parameters in Table I. The initial sample size was 50 or 10-times the dimensionality of the problem as it is a commonly suggested guideline to start with [34]. Simulations were then conducted to obtain the actual model responses for each initial sample point.

To generate meta-models, we used GP regression due to its effectiveness with limited datasets [35]. Each GP uses Bayesian principles and prior knowledge to create a posterior Gaussian distribution $f(x)∼GP( f ^(x), σ ^ 2(x))$⁠, where $f ^(x)$ is the predicted response (mean value) and $σ ^ 2(x)$ is the prediction variance. GPs utilize a kernel (or covariance function) to define the similarity between data points. Depending on the specific problem, various kernel functions can be selected or custom-designed for a GP (readers are referred to [36] for more detailed information on selecting a kernel/covariance function for a GP). The kernels used for the GP models are a combination of common choices, including radial basis function (⁠ $RBF$⁠), Matérn (⁠ $Mat$⁠), and dot product (⁠ $Dot$⁠) kernels, based on analyses explained in Appendix A 2.

For practical reasons, to perform parallel simulations, a batch of new samples was calculated according to the KB concepts [38]. The KB criterion updates the meta-models by believing (hallucinating) that the actual response of the new query is its GP mean-value prediction. This process of updating the quasi-Kriging model and applying the acquisition function is repeated until the desired number of samples within a batch (i.e., 10) is collected.

The simulations were run for each new batch and fed into the database to update the meta-models until a stopping criterion is met. To assess model performance and define a stopping criterion, the leave-one-out cross-validation (LOOCV) error can be expressed as a percentage error, from which the generalized mean squared cross-validation percentage error (GMSPE) is derived. Additionally, a test set of 27 data points was generated using Maximin distance design concepts. The actual responses of the test points were compared to those predicted by GP models, with percentage errors calculated and expressed as root mean squared percentage error (RMSPE). Figure 3 represents the GMSPE and RMSPE for (a) the BF response, the primary response of interest and (b) the MBL and MLL. As can be seen, both error criteria converge after a significant decline, indicating saturation and minimal likelihood for further improvement with additional iterations.

Conducting interpretability analysis is crucial for obtaining insight from developed machine-learning models. In this study, we applied SHAP interpretability analysis [25] to our developed meta-models to interpret the model predictions and explain how the system parameters (features) impact the model responses (predictions). The SHAP values are computed using a concept from cooperative game theory called Shapley values [39], which considers all possible feature combinations and determines the contribution of each feature to the model’s response. The SHAP analyses can assess feature importance models, depicting how features interact with one another and influence specific predictions. The interpretability analysis is discussed in further detail in Sec. III.

In this section, the SHAP analyses along with the results derived from the developed meta-models are discussed for the model responses, including BF, MLL, and MBL, using the final dataset of 120 data points.

Figure 4(a) illustrates the SHAP global feature importance for each system parameter or feature, where the absolute mean values are quantified for each individual system parameter. It is worth mentioning that the sum of all SHAP values for a given feature equals the difference between the actual prediction and the average predicted response. Therefore, larger SHAP values for a particular feature indicate a greater contribution to the deviation of the actual response from the average response, signifying its higher level of influence. The normalized gap (⁠ $H/ N k b k$⁠) emerged as the primary driver for the BF response, as evidenced by its larger corresponding SHAP value. This is also evident in the SHAP summary plot in Fig. 4(b), where the normalized gap is positioned at the top, distributing more widely. A more widespread spread of a feature’s points leads to a higher mean SHAP value, indicating a more significant impact on the model’s responses. Conversely, a narrow spread contributes to points clustered abound zero, reducing its mean SHAP value and showing little influence on the model’s responses. The gray scale intensity encodes feature values, with darker shades depicting higher feature values. The horizontal axis reflects the impact of these values on the BF response. Notably, higher values of the normalized gap (represented by black markers at the top row) are associated with lower predicted BFs (negative SHAP values). Consequently, as particle separation increases, the likelihood of chain bridging diminishes expectedly due to the reduced contact opportunities. Similarly, particle radius (⁠ $R p/ b k$⁠) influences the bridge formation as the second most important factor. Larger particles increase free space in the simulation box, leading to fewer bridges. This is further confirmed by the “effective cap fraction,” which measures the fraction of the surface area of the particle over which bridges are formed, as defined in our previous study [12]. Increasing particle radius at a constant gap and chain length reduces the effective cap fraction for bridge formation, leading simultaneously to increased availability of surfaces for loop formation. This phenomenon can also be observed in the Fig. 5, where the GP mean predicted bridge fraction is depicted against both $H/ N k b k$ and $R p/ b k$ for $ε s=8$  $k BT$⁠, $N k=10$⁠, and $N c/ N p=10$⁠, which we take as “standard” values. The other “standard value” is $R p/ b k=10$⁠, which values we assume in all plots except where otherwise noted. Time constants are always normalized by $τ 0= ζ b e a d b k 2/ k BT$⁠. Notably, BFs become negligible beyond a normalized gap of approximately 0.6, indicating a deficient bridge population and rare bridging dynamics. Under these conditions, most chains are expected to remain in loop configurations for highly extended periods, and any bridges that do form will break almost immediately, hence not contributing to any rheological effect of interest.

Figure 6(a) illustrates the relationship between the most important feature (⁠ $H/ N k b k$⁠) and the BF responses, revealing how the feature’s impact varies across different number of Kuhn steps (⁠ $N k$⁠). The plot is shaded according to $N k$⁠, with darker colors indicating higher $N k$ (larger chains) values and lighter colors representing lower $N k$ (shorter chains). Interestingly, as seen in Fig. 6(b), at normalized gaps below 0.2, longer chains correlate with increased BFs, whereas at normalized gaps greater than 0.2, longer chains exhibit lower BFs. As $N k$ grows, the rate of BF decline against $H/ N k b k$ becomes steeper. It should be noted that the normalized gap is the ratio of the closest distance between two adjacent particles divided by the chain length. Thus, to further demonstrate the dependency of BF on $N k$⁠, we calculated the available space for chains between particles. Defining the volume fraction (⁠ $ν$⁠) as the ratio of particle volume to total volume, the “available volume fraction” can be measured as $1−ν$⁠, accounting for both the normalized gap, $H/ N k b k$⁠, and $R p$⁠. Figure 6(c) depicts the predicted BFs vs both available volume fraction and $N k$ for systems of $ε s=8$  $k BT$⁠, $R p/ b k=10$⁠, and $N c/ N p=10$⁠. It is clear that larger chains are more likely to span larger gaps and, therefore, contribute to a higher BF at a constant $1−ν$⁠.

Conversely, the chain-to-particle ratio (⁠ $N c/ N p$⁠) and the attraction strength between stickers and colloidal particles (⁠ $ε s$⁠) exhibit minimal influence on BF model outcomes. Despite accounting for the excluded volume of beads, the relative abundance of chains and particles does not significantly affect BFs or other response metrics under equilibrium conditions. Thus, altering the chain population does not impact the proportion of chains forming bridges as shown in Fig. 7(a). Surprisingly, $ε s$ also demonstrates negligible impact on BF when other system parameters remain constant over the range of $ε s$ considered. Since $ε s$ governs the attraction between stickers and particles, increasing this attraction, for example, appears to equally stabilize both loop and bridge configurations—loops tend to retain their structure and bridges maintain their connections—leaving the bridges fraction nearly unchanged [as shown in Fig. 7(b)]. This is not surprising, at least for relatively high values of $ε s$ for which almost all stickers are bound to surfaces, which can be accomplished equally well through formation of either loops or bridges. It is the gap between particles that, therefore, controls the ratio of bridges to loops. However, it is evident that $ε s$ significantly influences the time a sticker remains bound to a particle, and hence $ε s$ strongly influences other system properties such as MLL and MBL, as will be discussed subsequently.

The SHAP global summary plots [Figs. 8(a) and 8(b)] reveal that the normalized gap is the most influential system parameter affecting the MLL response, as evidenced by its highest value in the SHAP global bar plot [Fig. 8(a)] and its wide distribution at the top of the SHAP summary plot [Fig. 8(b)]. The dark-colored dots on the right side of the SHAP summary plot demonstrate a positive correlation between the normalized gap and MLL, with increased gap values leading to a higher MLL. This suggests that there is an energetic preference for chains to remain in loop configurations rather than forming or transitioning to bridges, leading to an increased time spent in the loop configuration and consequently larger MLL values. As shown in Fig. 9, this phenomenon can be enhanced when the attraction strength (the second most influential parameter) also increases. Under the conditions of a large gap and strong attraction between stickers and particles, loops are further encouraged to maintain their configurations, leading to even larger MLL values.

Figure 10(a) illustrates the SHAP dependence plot of the normalized gap on MLL, considering its interaction with the number of Kuhn segments (⁠ $N k$⁠). When particles are positioned relatively close to each other, larger chains [represented by darker dots in Fig. 10(a)] do not necessarily exhibit a marked increase in MLL. This observation is also evident in Fig. 10(b), where MLL values are less sensitive to $N k$ at smaller normalized gaps than at higher gaps. This behavior is attributed to the increased likelihood of chains disrupting loops to form bridges at smaller gaps. Conversely, as the particle separation distance increases, the probability of transitioning from loop to bridge decreases, promoting more stable loops and consequently larger MLL values. This phenomenon is more pronounced for shorter chains.

Figure 8 suggests that the radius of particles (⁠ $R p$⁠) also positively contributes to MLL responses, where increased $R p$ leads to larger MLL values. As discussed previously, increasing $R p$ (while keeping other parameters constant) decreases the likelihood of the chains forming bridges, resulting in a higher proportion of chains maintaining loop configurations for more extended periods. An alternative perspective involves considering the available volume between the particles. Increasing the interparticle volume influenced by both the normalized gap and $R p$ is expected to correlate with larger MLL values. Figures 11(a) and (b) illustrate the relationship between predicted MLL values and available volume fractions for different $R p$ and $N k$ values, where clearly the higher available volume fractions correspond to longer MLL values. Interestingly, Fig. 11(b) suggests that at fixed available volume fractions shorter chains tend to exhibit longer MLL values due to the reduced probability of shorter loops traversing the same available space and transitioning to bridges compared to longer ones.

While the attraction strength between stickers and colloidal particles (⁠ $ε s$⁠) exhibited minimal influence on the BF in Figs. 4(a), 4(b), and 7(b), it emerges as a dominant factor affecting MBL as evidenced by its top ranking in Fig. 12(a). Figure 12(b) further emphasizes the importance of $ε s$⁠, with its positive contribution to MBL indicating that higher $ε s$ values contribute to longer-lived bridges and consequently larger MBL. Because forming bridges requires stretching chains to cross the gap, maintaining bridges costs stretching free energy, for which $ε s$ serves as a primary compensation force. After $ε s$⁠, $N k$ is the second most important factor influencing MBL. Longer chains can be stretched with less free energy cost and hence maintain their bridge configurations (once formed) longer than do shorter ones. As a consequence, higher or lower MBL values can be achieved by synergistically increasing or decreasing both $ε s$ and $N k$⁠, respectively. This insight is particularly useful for the design of coating fluid formulations.

Figure 13(a) illustrates the SHAP values for the normalized gap (⁠ $H/ N k b k$⁠) and its interaction with $N k$⁠, where the MBL values initially increases with increasing the $H/ N k b k$⁠, reaching a peak, followed by a decrease upon further increase in the gap ratios. The lower MBL values at smaller gaps can be attributed to the more confined space between particles, which pushes the stickers outside and places them in less optimal positions to maintain bridges. As the normalized gap widens, the influence of the normalized gap leads to decreased MBL values due to the increase energy cost of sustaining bridge configurations at larger separation distances. Moreover, at smaller normalized gaps, larger chains tend to increase MBL values. However, this correlation appears to be reversed as the normalized gap increases (darker colored dots correspond to lower SHAP values). This implies that the MBL values for chains with higher $N k$ decrease more sharply as the normalized gap grows. This behavior is further highlighted in Fig. 13(b), where the GP-predicted MBL is plotted against $H/ N k b k$ and $N k$ with other parameters at their standard values, as always. Notably, the position of the peak on the x-axis shifts rightward as the chain length increases. This trend can be attributed to the polymer end-to-end distance, which measures the distance between the two stickers within a chain. The MBL peak typically occurs at a distance between particles that approximates the equilibrium chain end-to-end distance, which is approximately $H/ b k N k$⁠, representing the most favorable distance for sticker interaction under equilibrium as confirmed by the peak in MBL at $H/ b k N k$ around unity in Fig. 13(c) for different $N k$⁠. This behavior provides valuable insight that can guide the development of design principles for formulating waterborne coating fluids.

This study explored the association dynamics of telechelic polymers with colloidal latex particles, a pivotal factor in understanding and hence controlling the rheology of waterborne coating fluids. Recognizing the limitations of existing simulation techniques in capturing the full spectrum of relevant time and length scales, we developed a novel active learning-based approach to bridge the gap between high-resolution models and those with lower resolution. Initially, we developed BD simulations as precise yet computationally expensive models, resolving the fine-grained interaction of telechelic chains with colloidal particle surfaces. As a result, we could measure the system’s microscopic responses, such as the bridges and loops lifetimes and the BFs. These simulations allowed the generation of initial and sequential system responses for any points within the design space, utilizing an active learning approach. The active learning meta-models developed here are capable of predicting the lifespans of bridges and loops and the populations of bridges within the design space. This achievement is realized with a notably reduced number of samples required to map the system parameters to the corresponding responses from thousands to one hundred samples. Moreover, the interpretability of these meta-models enhances the model transparency and further facilitates refinement of formulations by focusing on the most influential design parameters.

Analyzing the meta-models provides insights into the interplay between system parameters in predicting MLLs and MBLs, and BFs. By adjusting parameters such as polymer length (represented in Kuhn steps), the number of chains per particle (polymer concentration in the formulation), particle spacing (which reflects the particle volume fraction in the solution), and the attraction strength between stickers and particles (number of the methylene in the hydrophobic end caps), experimentalists can tailor formulations to fine-tune the rheological behavior of waterborne coatings. The findings indicate that the normalized gap between particles plays a prominent role in bridge formation and loop lifetime. Conversely, the attraction strength between the stickers and particles has a minimal impact on the BF. Furthermore, this study identified a synergistic effect of chain length and attraction strength on promoting extended bridge lifetimes, resulting in longer relaxation times. This has been found experimentally [40,41] to increase the viscosity and stability of analogous systems, i.e., interconnected micelles by bridging chains. The developed meta-model effectively demonstrates that the optimal distance between particles maintaining bridges over an extended period corresponds to the end-to-end distance of the chains under equilibrium conditions, which refers to a more stable coating suspension in storage. In addition, the intelligent search strategy developed in this work can easily be adopted to accelerate and guide the design of laboratory experiments and is not limited to computational studies.

This study has laid a solid foundation for future research on the interaction between telechelic polymers and colloidal latex particles. Several avenues can be explored to extend and deepen our understanding of this complex systems. Integrating these meta-models with coarser-grained simulation methods, such as Pop-BD [12], could provide a more comprehensive understanding of the system at different lengths and time scales. The significant advantage of this research lies in its ability to formulate quantitative models suitable for integration into Pop-BD simulations. In addition to the gap-dependent BF, Pop-BD simulations require the loop-to-bridge transition rate (L) and the bridge-to-loop transition rate (M) that are equivalent to 1/MLL and 1/MBL, respectively. These Pop-BD input values are represented for the standard values (i.e., $ε s=8$  $k BT$⁠, $N k=10$⁠, $R p/ b k=10$⁠, and $N c/ N p=10$⁠) in Fig. 14. By providing accurate estimates of BF and these transition rates, our predictive meta-models allow the Pop-BD simulations to accurately capture the higher relaxation times of systems with larger colloidal particles, which is essential for bridging the microscale interactions to macroscale rheological behavior. With the developed meta-modes, our future research will focus on enhancing the predictive power of Pop-BD simulations informed by developed meta-models and gaining a more comprehensive understanding of the rheological properties at the macroscale.

Future work should also focus on extending the meta-models to larger values of the particle size, chain length, and numbers of chains per particle, and sticker strength. This could be done by a combination of additional runs and using ML methods to extrapolate from smaller parameter values to larger, using more sparse sampling for the more expensive simulations to conserve computer power. Moreover, investigating the system’s behavior under shear conditions would provide valuable insights into the rheological properties of waterborne coatings. Refined gap-dependent analyses of these systems under shear forces are of great interest, as highlighted in [42]. These future studies ultimately pave the way for optimizing waterborne coating formulations to achieve the desired macrolevel rheological properties.

Consequently, this work represents a significant step forward in the multiscale modeling of waterborne coatings, offering a computationally efficient method to provide detailed microscopic-level information for use in macroscopic simulations. Offering a comprehensive understanding of how various factors such as polymer chain length, particle size, and attraction strength contribute to association dynamics under equilibrium conditions, these meta-models have minimized the need for extensive trial-and-error simulations and can accelerate the development of more efficient and environmentally friendly coatings.

This research was supported by The University of Melbourne’s Research Computing Services and the Petascale Campus Initiative. J.A. acknowledges the training facilities at the Australian Research Council (ARC) Industrial Transformation Training Centre in Optimisation Technologies, Integrated Methodologies, and Applications (OPTIMA). R.G.L. is grateful for the financial support of the National Science Foundation of the USA under Grant No. DMR-2100513. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NSF.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The developed BD simulations update the positions of all beads and stickers in the system at each time step, enabling the study of their dynamics over time. A time step of $5× 10 − 5[t]$ was used, consistent with the literature [27]. As time progresses, the positions of the stickers are tracked at specific time intervals, known as dumped time steps. These time steps are typically set to be as low as possible to capture accurate particle positions. To manage the volume of data generated, a dumped time step of 10 was used, ensuring that it is significantly smaller than any bridge or loop lifetime, as calculated according to the methods described in Sec. II A. It is crucial to allow simulations to run for a sufficiently long duration to capture adequate statistics of bridge and loop lifetimes before calculating averages. Failing to do so or imposing a cap on the total run time can lead to inaccurate results, particularly for loop lifetimes, which can be significantly longer than the capped total run time. Given the random and independent nature of loop and bridge lifetimes, their observed distributions should closely resemble a Poisson distribution [43]. The histogram plots of the samples are presented in Fig. 15 for a system with $H/ N k b k=0.2$  $N k=10$ of $ε s=8$  $k BT$⁠, $R p/ b k=10$⁠, and $N c/ N p=10$⁠. As observed, the frequency of occurrence of loop or bridge samples decreases exponentially as the duration of events increases.

We used the criteria mentioned above to assess the kernel performance and select the best kernel for generating new samples. After generating the initial samples, the LML values and GMSPE were calculated for different kernel combinations at each iteration. These results were used to compare and evaluate the kernels that have the maximum value for LML and the minimum value for GMSPE. Another criterion considered was the simplicity of the kernels or having fewer hyperparameters to avoid overfitting. Once a satisfactory level of accuracy was achieved (below 5 $%$ in GMSPE computed for the BF responses), the final selected kernels were retained for the remaining iterations. Figure 16 represents the GMSPE and LML values for different combinations of kernels used to predict the three responses: BF, MLL, and MBL. The data are based on a population of 80 samples. According Fig. 16(a), $( M a t r n ( ) + D o t P r o d u c t ( ) 2 ) 2$ was selected because it has a relatively low GMSPE value (although not the lowest), a high LML value, and the fewest number of hyper parameters. Similarly, $(Matrn()×RBF()$ $× D o t P r o d u c t())$ and $( M a t r n ( ) + D o t P r o d u c t ( ) 2 ) 2$ were selected to represent kernels for MLL and MBL, respectively.