Open Access Article
Lijie
Ding
a,
Chi-Huan
Tung
a,
Jan-Michael Y.
Carrillo
b,
Wei-Ren
Chen
a and
Changwoo
Do
*a
aNeutron Scattering Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA. E-mail: doc1@ornl.gov
bCenter for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
First published on 23rd June 2025
We develop Monte Carlo simulations for uniformly charged polymers and a machine learning algorithm to interpret the intra-polymer structure factor of the charged polymer system, which can be obtained from small-angle scattering experiments. The polymer is modeled as a chain of fixed-length bonds, where the connected bonds are subject to bending energy, and there is also a screened Coulomb potential for charge interaction between all joints. The bending energy is determined by the intrinsic bending stiffness, and the charge interaction depends on the interaction strength and screening length. All three contribute to the stiffness of the polymer chain and lead to longer and larger polymer conformations. The screening length also introduces a second length scale for the polymer besides the bending persistence length. To obtain the inverse mapping from the structure factor to these polymer conformation and energy-related parameters, we generate a large data set of structure factors by running simulations for a wide range of polymer energy parameters. We use principal component analysis to investigate the intra-polymer structure factors and determine the feasibility of the inversion using the nearest neighbor distance. We employ Gaussian process regression to achieve the inverse mapping and extract the characteristic parameters of polymers from the structure factor with low relative error.
To understand the structure and behavior of the charged polymers, both experimental and theoretical approaches have been employed. Experimental techniques such as small-angle scattering12 (SAS) including X-ray scattering13 and neutron scattering14,15 have proven indispensable for understanding these properties of the charged polymers.16 Scattering methods provide insights into the nanoscale structure and dynamics of charged polymers, enabling the characterization of key conformational parameters such as radius of gyration, persistence length, and inter- and intra-molecular interactions. Theoretical and computational approaches, including analytical models17,18 and computer simulations, complement experimental efforts by capturing the fundamental physics of charged polymer systems. Techniques such as molecular dynamics19,20 (MD) and Monte Carlo21,22 (MC) simulations have provided significant insights into polymer configurations, bending rigidity, and electrostatic interactions.
Despite the progress made on both the experimental and theoretical fronts, bridging the scattering function measured in SAS experiments with the polymer parameters used for modeling charged polymers in theory and simulations remains a significant challenge. The difficulties lie in extracting physical quantities about polymer conformation by decoding the scattering function. Recent advances in machine learning (ML) have opened new avenues in scattering analysis, enabling parameter extraction without requiring explicit analytical forms of the scattering function.23 By training ML models on simulation-generated data, it becomes possible to establish an inverse mapping from the scattering function to the underlying model parameters. This approach has shown promise in a variety of systems, including colloids,23–25 polymers,26–29 and lamellar structures.30,31 These applications demonstrate the potential of ML to bridge the gap between experimental scattering data and theoretical models, providing a robust framework for parameter extraction in complex systems.
In this work, we introduce such an inversion by the ML approach for the charged polymer system, where the data are generated using MC simulations. The polymer configuration is governed by the intrinsic bending stiffness, charge density and salt concentration of the surrounding medium. We first investigate the effects of these key variables on polymer conformation and then calculate the intra-polymer structure factor. To assess the feasibility of inversion, we perform principal component analysis of the scattering data and quantify the feasibility using the nearest neighbor distance of the polymer parameters in the structure factor space. Finally, we employ Gaussian process regression (GPR) to extract both the conformational and energy-related parameters of the polymers from the structure factor, demonstrating the accuracy and robustness of this approach.
![]() | (1) |
is the Yukawa potential, or screened Coulomb potential,20,32 that models the charge interaction, A is the interaction strength between charged monomers, λD is the Debye screening length,33 and
is the distance between joints i and j. In addition, the self-avoidance of the polymer is enforced by adding hard sphere interaction of diameter lb between different joints. The interaction strength
is directly related to the charge density of the polymer σe, where ε is the dielectric constant of the medium. The Debye screen length
, where kB is the Boltzmann constant, T is the system temperature, e is elementary charge, and
is the ionic strength, in which zi and ni are the charge number of the number density of ion species i, respectively.
is updated using two MC moves: continuous crankshaft and pivot. Crankshaft picks two random joints on the polymer chain and rotates all the bonds between them for a random angle within the interval [−ϕc, ϕc]. Pivot randomly selects one joint k on the chain and rotate the preceding sub-chain (k, …, N) within a cone of angle ϕp(k) centering at the original orientation. To improve the acceptance rate of these updates and thus boost the efficiency of the simulation, the crankshaft rotation angle is adjusted according to the bending modulus such that
, and the pivot rotation angle ϕp = ϕc. Combining these two moves allows full exploration of the polymer configuration with the contour length fixed and the polymer conformation calculated using this algorithm has been benchmarked against theoretical calculations. More details on the MCMC simulation can be found in our previous paper.34
To better characterize and understand the conformation of the charged polymer, we calculate the radius of gyration, bond angle correlation and structure factor of the polymer. The radius of gyration square is
, where the 〈…〉ij denotes the average of all pair of joints. The bond–bond correlation is 〈cos(θ(s))〉 = 〈
i·
i+s〉i where 〈…〉i denotes the average over all bonds and s represents the contour distance between two bonds along the polymer chain. Finally, the isotropic intra-polymer structure factor12,14 is given by:
![]() | (2) |
S(q) and carry out principal component analysis for the data sets. The S(q) is calculated for 100 q ∈ [10−1, 1], uniformly placed in the log scale, and κ ∼ U(5, 50), A ∼ U(0, 10) and λD ∼ Ud(1, 10), where U(a, b) is the uniform distribution in the interval [a, b] and Ud(a, b) is the discrete uniform distribution. Similar to previous work,23 we use singular value decomposition (SVD) to find the three most important bases of the 4000 × 100 matrix F = {log
S(q)}, such that F = UΣVT. The diagonal entries of Σ2 are proportional to the weight of the variance of the projection of F onto each principal vector of V. Projecting F onto the first few bases provides a way to analyze F in a dimensionally reduced space. A useful tool to study the distribution of the polymer parameters Y = {(κ,A,Rg2/L2,R2/L2)} is to calculate the nearest-neighbor distance of ζ ∈ {κ,A,Rg2/L2,R2/L2} on the F manifold. For n-number of vectors, x1, x2, …, xn, the first nearest neighbor is defined as NN1(xi) = argminxj≠xi|xj − xi|; similarly, the second nearest neighbor is NN1(xi) = argminxj≠xi,NN1(xi)|xj − xi|, and we define the normalized nearest neighbor distance DNN for the ζ(x) as:![]() | (3) |
S(q) can be traced back to unique changes in ζ. Concretely, large DNN(ζ) indicates that minor variances in log
S(q) can map to large jumps in ζ, signaling regions where inversion is unstable or degenerate. Whereas small DNN(ζ) means that even significant noise in log
S(q) produces only modest shifts in ζ, thus the inverse mapping remains well-conditioned and robust.
S(q), to the system parameters, or inversion targets y = (κ, A, Rg/L2, R2/L2), we employ a Gaussian Process Regression (GPR) model trained on data generated through Monte Carlo (MC) simulations. Under the framework of GPR,36,37 the goal is to obtain the posterior distribution p(Y*|X*, X, Y) for the function output y. In this setup, the training and test sets are defined as X = {log
S(q)}train and X* = log
S(q)test, respectively, while Y and Y* correspond to the inversion targets (κ, A, Rg/L2, R2/L2). GPR assumes a Gaussian process prior over the regression function, g(x) ∼ GP(m(x), k(x, x′)), where m(x) is the prior mean function and k(x, x′) is the covariance kernel. The joint distribution for the Gaussian process is expressed as follows:![]() | (4) |
Here, we use a constant prior mean function m(x), while the kernel function is modeled as a combination of a Radial Basis Function (RBF) and a white noise term:
![]() | (5) |
S(q)}. We then discuss the feasibility of inversion based on the SVD of F. With the feasibility established, we finally test our trained GPR for the inversion.
When the polymer is only subjected to bending κ, or in the case of A = 0, the polymer is a classic semiflexible polymer whose bond angle correlation can be described by a single exponential decay:
〈cos θ(s)〉 = e−s/λ0 | (6) |
〈cos θ(s)〉 = (1 − α)e−s/λ1 + αe−s/λ2 | (7) |
Fig. 2(a) shows the bond angle correlation function 〈cos
θ(s)〉 for various screening length λD, and the fitted lines are calculated using the single scale model as in eqn (6). As λD increases, the single scale model fitting starts to diverge from the data point, indicating the necessity of switching to the double length scale model (eqn (7)); Fig. 2(b) shows such fitting results, and the two length scale model can still describe the decay of 〈cos
θ(s)〉 at large λD.
![]() | ||
Fig. 2 Different length scales of the charged polymer, fitted using both single length scale and double length scale models. (a) Bond angle correlation 〈cos θ(s)〉 for various screening length λD with κ = 30, A = 5, solid lines are fitted using a single length scale (eqn (6)). (b) Similarly, but fitted using double length scale (eqn (7)). (c) Three persistent lengths, λ0 for the solid line, λ1 for the dashed line and λe for the dotted line, versus screening length λD for various κ with A = 5. (d) Similar to (c), but for various A with κ = 30. | ||
Fig. 2(c) show all three length scales λ0, λ1 and λeversus screening length λD for various bending stiffness κ. At low λD the one length scale still fits the bond angle correlation data, and increases with increasing λD. When switching to the two length scale model, the long length scale λ1 increases with increasing λD, while the short length scale λe decreases and deviates from λ1 and then plateaus. The plateau value increases with bending stiffness κ. The switch from the one-length scale λ0 to two-length scale (λ1, λe) in the plot is determined by monitoring the divergence between the λ0 and λ1 when fitting the correlation function at low screening length λD. Fig. 2(d) shows a similar result but for various charge interaction strength A. Similar to its effect on the end-to-end distance and radius of gyration, A amplifies the effect of increasing λD, while the short length scale λe plateaus at a similar value for various A, confirming it corresponds to the bending stiffness κ.
, with all bonds pointing to the same direction. Fig. 3(a) shows the variation of structure factor S(q) for various bending stiffness κ. Compared to the solid rod, the polymer structure factor shows a bump at a structure vector q range comparable to its radius of gyration. Fig. 3(b) shows the structure factor of the polymer divided by the rod S(q)/Srod(q), where the bump is better shown. As the bending stiffness κ increases, the peak in S(q)/Srod(q) lowers and the corresponding q value also decreases, indicating an increase of the characteristic length. Fig. 3(c) and (d) shows the S(q)/Srod(q) for various charge interaction strength A and screening length λD, and both show similar effects on the structure factor of the polymer as they make the polymer more extended and stiff.
To better analyze the structure factor of the charged polymer, we carry out principal component analysis described in Sec. 2.3. By decomposing the F = {log
S(q)} into F = UΣVT, we find that the singular value Σ decays rapidly versus its rank, as shown in Fig. 4(a), indicating we can represent the log
S(q) ∈ F using few bases. Fig. 4(b) shows the first 3 singular vectors, and Fig. 4(c) shows the projection of a structure factor S(q) onto each basis, and the reconstruction from only the 3 bases closely matches the original S(q). This decomposition will allow us to further determine the feasibility of extracting these polymer parameters from the structure factor.
S(q) ∈ F into the space spanned by the first 3 singular vectors (V0, V1, V2), and the corresponding 3 coefficients of each log
S(q) correspond to a single point in the
space. As shown in Fig. 5(a–c), the end-to-end distance R2/L2, radius of gyration Rg2/L2 and bending stiffness κ are all well spread out on in the FV manifold, indicating they are eligible to be extracted from the structure factor. Fig. 5(d) shows the distribution of charge interaction strength A and it is unclear if it can be extracted due to some randomness in the distribution.
Intuitively, when then screening length λD is very small, the effect of the charge interaction becomes negligible, preventing it from having a meaningful impact on the structure factor S(q), thus it is not expected to have A feasible for extraction from the S(q) at low λD. To quantify this feasibility, we slice the structure factor data set F = {log
S(q)} into different slices for different screening lengths λD, and calculate the nearest neighbor distance for each slice. As shown in Fig. 6(a), we plot 3 slices of the charge interaction strength A distribution, and the randomness reduces as the screening length λD increases. Quantitatively, Fig. 6(b) shows the nearest neighbor distance DNN for each polymer parameter and DNN(A) is much larger than that of the others when the screening length λD is small, and then it decays to lower value as the λD increases, leading to a more significant impact of the charge interaction strength A on the polymer conformation. This indicates the charge interaction strength A, which is directly related to the charge density of the polymer, is still extractable if the screening length is large enough.
![]() | ||
| Fig. 6 Nearest neighbor distance analysis of the charge interaction strength A. (a) Value distribution of A in the SVD space for various slices of screening length λD, the axes are the same as in Fig. 5. (b) Nearest neighbor distance DNN for various polymer parameters versus different slices of the data F separated by the λD value. | ||
S(q)}as the training set {log
S(q)}train, and then test the trained GPR using the remaining 30% data {log
S(q)}test by comparing the actual polymer parameters with the ones extracted from the structure factor S(q). The split between the training and testing data is random. To obtain the trained regressor, we need to find the optimized hyperparameters (l, σ) for each inversion target, or polymer parameters. We search for the (l, σ) that maximize the log marginal likelihood,36 which are shown in Fig. 7.
Fig. 8 shows a comparison between polymer parameters ((R2/L2,Rg2/L2,κ,A)) obtained from ML inversion and the corresponding reference used in or calculated through MC simulation. We note that due to the high nearest neighbor distance DNN(A) of charge interaction strength at low screening length λD, we only used data with λD ≥ 4 for the inversion of A. Nevertheless, the data agree well, and lie closely along the diagonal line, with relatively low error, and for polymer parameter ζ, the relative error between MC reference ζMC and ML inversion ζML is estimated by Err = 〈|ζMC − ζML|/max(ζMC,ζML)〉, where 〈…〉 is the average over all data points. The relative error is annotated on each panel of Fig. 8 and shows very high precision for ((R2/L2,Rg2/L2,κ) and good precision for A. While the errors for end-to-end distance R2, radius of gyration Rg2 and bending modulus κ are very small, the error for charge density is relatively large as we are including data with all screening length λ ≥ 4. In practice, the screening length can be estimated based on the solvent conditions; a reduced range of λD will lead to better accuracy in the extraction of charge density A.
Our approach provides a unique method to obtain the bending stiffness and the charge density σe, which is directly related to the charge interaction strength
using the scattering data. A natural next step would be to carry out a SANS experiment for some charged polymer sample, and apply our approach on the experimentally measured SANS data. In practice, this approach assumes the experimental data falls within the range of training data, and a procedure of trial and error maybe required based on the fitting results, in which the training set needs to be expanded as needed. In addition, experimental data naturally come with noise, for which a denoising procedure40 can be helpful, and the analysis of noisy data will naturally provide uncertainties by the GPR.28 Moreover, the effect of noise for the GPR prediction can be systematically studied by measuring the accuracy of the inversion when different levels of noise are added to the testing data. Finally, this framework can be expanded to the study of more complicated charged polymer systems including charge-patterned polypeptides,41 alternating copolymers42 and zwitterionic patterned polymers.43 To study these systems, it is required to model the polymer energy accordingly. It is natural to introduce variable charge interaction strength A for different monomer segments to model the charge pattern and polarity, and a screened dipole–dipole interaction can be used for modeling the zwitterionic polymer.
| This journal is © The Royal Society of Chemistry 2025 |