Short-time dynamics of lysozyme solutions with competing short-range attraction and long-range repulsion: Experiment and theory

Experimental and theoretical studies of sophisticated bio-particle systems have become an important part of soft matter science. In particular, the understanding of concentrated protein solutions has attracted much attention due to its importance for both biological systems and pharmaceutical industrial products.^1–17 However, arguably, it is still extremely challenging in colloidal science to accurately describe the structure and dynamics of concentrated protein solutions, as many proteins can have irregular shapes, complex protein-protein interactions, and heterogeneous surface properties. Therefore, much effort has been devoted to study globular proteins which can be approximated as spherical objects with isotropic interaction potentials,^2,5,6,9,12 through which useful insights can be obtained for more complex protein systems, such as monoclonal antibodies.^14–16 Among the various properties of protein solutions, there is recent interest in understanding the short-time diffusion coefficients, d_s, and the viscosity, η, of concentrated protein solutions for both academic research and industrial applications.^{4–6,8,11,13} The short-time diffusion coefficients obtained by neutron spin echo (NSE) and backscattering techniques have been useful for understanding the cluster formation of protein solutions and corresponding viscosity.^4,7,13,17 The prediction of viscosity of complicated protein solutions is critically important for the pharmaceutical industry to develop more effective drug delivery strategies for some cancer treatment therapeutic proteins.^14,16,17

Calculating d_s and η has been extensively studied for colloidal systems with hard-sphere interaction. The viscosity behavior of hard-sphere systems can be satisfactorily described by mode coupling theory (MCT).^18–21 The change of d_s in hard-sphere systems has been investigated both experimentally²² and theoretically.²⁰ Short-time dynamics of charged colloidal systems with electrostatic repulsions have been examined using existing theories, which is compared with accelerated Stokesian dynamics (ASD) simulations to investigate the applicability and accuracy of various theoretical expressions.²³ 

However, despite the importance in calculating d_s and η for protein solutions, there are very limited studies of using colloidal theories to understand protein solutions. At low salt concentrations, it has been shown that the interaction potential between bovine serum albumin (BSA) proteins can be well described by a hard-sphere potential with an additional electrostatic repulsion.⁹ At this condition, the short-time diffusion coefficient of BSA proteins has been studied by quasi-elastic neutron scattering and is shown to follow the prediction of a system with purely hard-sphere interaction.¹³ The experimental viscosity of BSA protein solutions has been used to examine different theories.⁸ 

However, for many proteins in solution, it is common that the interaction contains both a short-range attraction and long-range electrostatic repulsion (SALR), such as in lysozyme and some monoclonal antibody protein solutions at low salt concentration.^1,5,11,14 SALR systems demonstrate interesting solution structures, as the competition of the attraction and repulsion can introduce intermediate range order (IRO).^5,24 Depending on different combinations of the attraction and repulsion contributions in protein solutions, different types of IRO structures can be introduced.^11,25–28 When the attraction has a very short range, which is the case for many protein systems such as lysozyme, a recent work has shown that there can be different liquid states including dispersed fluid states, clustered fluid states, random percolated states, and cluster percolated states.²⁹ 

Despite many interesting studies on SALR systems, there are very few studies on the short-time diffusion coefficients and viscosity behavior. In fact, studies on SALR systems of model colloidal particles have used micrometer sized particles and have mostly focused on the solution structures and the gelation/glass transitions.^25,26 In contrast, the experimental-theoretical study of protein solutions with a SALR interaction has mostly focused on short-time diffusion coefficients and solution viscosities,^{4–7,12,14,16} with the salient hydrodynamic interactions (HIs) between the proteins being disregarded or strongly approximated in the theoretical treatments.

Only very recently, some of the authors of this paper have theoretically investigated the short-time dynamics of SALR systems with two-Yukawa potential interactions and HIs included.³⁰ This theoretical work for the first time demonstrated that an intermediate range order (IRO) peak can appear in the hydrodynamic function. Moreover, an unexpected non-monotonic temperature dependence of the mean particle sedimentation velocity was predicted theoretically for homogeneous systems.

In this paper, we test the previously developed theory to study the short-time dynamics and viscosity of lysozyme solutions, where the interaction has to be described by a short-range attraction and long-range repulsion. Using the experimentally obtained interaction potential parameters, we have first identified that the studied lysozyme samples are mostly in a dispersed fluid phase state, while some of the samples at low temperature and high concentrations are in the randomly percolated state. Using the corresponding structure factors, the hydrodynamic function, H(q), short-time diffusion coefficient D(q), and viscosity, η, are calculated and compared against experimental values for lysozyme solutions in the dispersed fluid state. We critically examine the accuracy and limitations of the current theoretical model applied to lysozyme samples. This work also provides theoretical insights into the peculiar experimental features typically observed in SALR systems.

The theoretical methods³⁰ and experiments^11,31 employed in this work have been reported in publications separately by some of us. Brief overviews are provided here for clarity and completeness since in this paper the theory is compared against the experimental results for the first time.

Model interaction potentials such as the two-Yukawa^{3,5,28–30,32–34} and the generalized Lennard-Jones Yukawa (LJY)^6,29,35,36 potentials, respectively, have been widely used to describe the inter-particle potential in different phase states of SALR systems including the dispersed fluid,^{4,6,7,29–31,36} clustered-fluid,^{4,6,7,29,31,32,36,37} random percolated,^6,29,31,32 and glassy states.^{6,11,25,26,31} In the present work, we use the hard-core plus two-Yukawa (HCDY) pair potential to represent the interaction for a SALR model system. The HCDY potential can be expressed as

where x = r/σ is the inter-particle center-to-center distance, r, in units of the particle diameter σ, and β = 1/k_BT. Moreover, z₁ and z₂ determine the range of the attractive and repulsive Yukawa potential parts in units of σ, respectively, and K₁ and K₂ are the respective short-range attractive and long-range repulsive potential strengths in units of k_BT. The depth of the attractive well (or the net attraction) is given by βV (x = 1⁺) = K₂ − K₁.

Microstructural properties of lysozyme samples are reproduced using Monte Carlo (MC) computer simulations^38,39 consisting of 1728 particles in the NVT ensemble within the one-phase region and by employing periodic boundary conditions. Proteins are represented as spherical particles interacting isotropically via the hard-core plus two-Yukawa (HCDY) potential. Starting from a simple cubic lattice, each system is thermally equilibrated for 2 × 10⁷ steps, after which thermodynamic and structural parameters were averaged over 4 × 10⁴ independent configurations. The initial displacement distance of 0.1, where all distances are normalized by the particle diameter σ, is dynamically adjusted to maintain an acceptance ratio of 30%. The observables were averaged over 10 different seeds to reduce the intrinsic uncertainties. In all cases, a system size of N = 1728 particles was sufficient to avoid artificial size effects and optimize the computational time.

The microstructure is represented by first calculating the radial distribution function, g(r), and static structure factor, S(q), for each considered system. The radial distribution function is calculated by averaging all particle configurations using direct summation according to^38,39

where δ is the three-dimensional delta function, N is the number of particles in the system volume V, r is the distance vector of modulus r, and r_ij is the vector distance between the centers of particles i and j. Furthermore, ⟨⋯⟩_eq denotes an equilibrium average. The protein solution structure factor is then determined by its relationship with g(r),

where ρ = N/V is the particle concentration, q = (4π/λ)sin(θ/2) is the modulus of the scattering wave vector, and λ and θ are the incident wavelength and scattering angle, respectively. Cluster formation is identified when a particle is less than the cut-off distance, r_c, away from a neighboring particle, which is defined as the separation where the interaction potential produces a zero interaction energy.²⁹ Connectivity calculations yield an account of all cluster sizes, which is summarized in a cluster number distribution, n(s). The cluster number distribution is normalized by the cluster size, s, and system size, N_p, resulting in the cluster size distribution (CSD),^29,32

where N_p is the total number of particles in the simulation box such that $∑s=1NpN(s)=1$⁠. N(s) represents the fraction of particles contained in clusters of size s and is a normalized function for all cluster sizes. This normalized CSD defines the state of the fluid at each set of conditions using definitions established previously.²⁹ Along with the CSD, a nearest neighbor number distribution (NND) function, f(N), is calculated as the average fraction of particles having N neighboring particles whose centers are located inside the attractive well range of the HCDY potential.

For the calculation of equilibrium pair functions of SALR systems using liquid-state integral equation theory methods, the importance of using a self-consistent closure relation was shown recently in Ref. 40. A self-consistent hybrid scheme for systems with attractive and repulsive pair potential parts is used here based on the method proposed by Zerah and Hansen (ZH).⁴¹ The ZH scheme interpolates between the hypernetted chain closure (HNC) for long and the soft-core mean spherical approximation (SMSA) for short particle pair distances r.

The ZH scheme is particularly well suited for systems with soft-core repulsion and an attractive interaction part.^41–43 To describe such SALR systems theoretically, we split the total pair potential V(r) into a reference part, V₁(r), and a perturbation part, V₂(r), selected here as

and

respectively, where r_shift = r_max > σ, and where r_max is the pair distance at which V(r) is at its maximum. Accordingly, V₂(r) is purely attractive, while V₁(r) is purely repulsive (c.f. Ref. 44). The ZH closure reads

with the mixing function

Here, the mixing parameter ϑ is determined self-consistently from enforcing equality of the isothermal compressibilities derived using the virial pressure and compressibility routes, respectively. In taking the concentration derivative of the virial pressure, we assume for simplicity that the mixing parameter is density independent. This approximation is justified since ϑ is only a weakly varying function of the volume fraction ϕ.⁴¹ In the limit of ϑ → ∞, the ZH closure reduces to the hypernetted chain closure (HNC), while in the opposite limit ϑ → 0, the soft MSA closure relation given by

is recovered.

Since there exist different choices for V₁(r) and V₂(r), the comparison with simulation data for the radial distribution function, g(r), is a necessary prerequisite to assess the accuracy of the ZH scheme.⁴¹ For the dispersed-fluid phase, the ZH approach with this splitting of V(r) yields results for g(r) in excellent agreement with computer simulation predictions based on the MC, MD and multi-particle collision (MPC) simulation methods, as it is shown in Refs. 30, 36, and 44.

Short-time diffusion in colloidal suspensions is commonly assessed experimentally by measuring the intermediate scattering function S(q, t), where

is an exponentially decaying function for correlation times, t, small compared to the structural relaxation time τ_d = a²/d₀. Here, d₀ is the Stokes-Einstein single-particle translational diffusion coefficient for a spherical particle (protein) of radius a = σ/2 in a solvent of viscosity η₀, given by d₀ = (k_BT) /(6πη₀a) for stick surface boundary conditions. For t ≪ τ_d, the configuration of particles is hardly changed by diffusion such that the particle dynamics is influenced solely by hydrodynamic interactions (HIs). In Eq. (10), D(q) = d₀H(q)/S(q) denotes the short-time diffusion function, and H(q) is the positive-valued hydrodynamic function. The latter characterizes the influence of HIs on short-time diffusion, and it can be calculated in overdamped Brownian dynamics starting from

The thermodynamic limit, denoted by $lim∞$⁠, is taken here to describe a macroscopically large scattering volume. The μ_ij(X) are the protein configuration (i.e., $X=r1,…,rN$⁠) dependent mobility matrix tensor elements linearly relating the hydrodynamic force on a particle j to the velocity change of a particle i owing to the solvent-mediated HIs. Moreover, k_BT μ₀ = d₀ and $q^=q/q$⁠. The diagonal terms in Eq. (11) for which i = j give the wavenumber-independent short-time self-diffusion coefficient d_s, while the off-diagonal (i ≠ j) terms sum up to the wavenumber-dependent distinct hydrodynamic function part, H^d(q), of H(q), that characterizes the hydrodynamic force-velocity cross correlations.

For the calculation of H(q), we use the so-called hybrid Beenakker-Mazur pairwise additivity (BM-PA) scheme which combines the well-established Beenakker-Mazur (BM) effective medium method,^45,46 where many-body HIs are approximately included for the calculation of H^d(q), with the Pairwise Additivity (PA) approximation of HIs used for the calculation of the self-part, d_s, of H(q). In recent work by some of the present authors,³⁶ the good accuracy of the BM-PA scheme for the calculation of short-time diffusion properties of SALR systems has been shown by the comparison with elaborate multi-particle collision (MPC) dynamics simulation results where HIs are fully accounted for. For details about the employed BM-PA hybrid method, we refer to Refs. 30, 36, 44, and 47. We emphasize that for a given set of HCDY pair potential parameters, the BM-PA theory is predictive and requires no fitting to the data. The HCDY potential parameters are determined by fitting the small angle neutron scattering (SANS) data using the ZH integral equation scheme in conjunction with the decoupling approximation.

In addition to short-time diffusion properties, we have determined theoretically the zero-frequency low-shear viscosity η = η_∞ + Δη as a long-time transport property. The high-frequency viscosity part, η_∞, is of purely hydrodynamic origin and quite accurate analytic tools for its calculation are available. In contrast, the calculation of the shear relaxation part, Δη, is more demanding because one needs to account for the shear-induced deformation of next-neighbor cages influenced by direct and hydrodynamic interactions alike. To calculate Δη, we use a simplified mode-coupling theory (MCT) expression where the viscosity is obtained in a first iteration step of the self-consistent MCT equations, by relating η to the time evolution of S(q, t). Explicitly, we use the expression^20,47

where y = qσ and S′(y) = dS(y)/dy. In this expression, the contribution of the HIs to the MCT shear relaxation vertex function is omitted, which can be partially justified by the fact that the associated hydrodynamic mobility tensors relating shear strain to hydrodynamic particle force dipoles (stresslets) are rather short-ranged (see, e.g., Ref. 21). HIs enter in Eq. (12) only through S(q, t), approximated by its short-time form

valid for t ≪ τ_d. Consequently, the relaxation part Δη is underestimated by the first-order MCT as compared to the fully self-consistent MCT. This underestimation is expected to become more pronounced at larger particles volume fraction ϕ = (π/6)ρσ³.

For calculating the high-frequency viscosity contribution, η_∞, we use the PA method where two-body HI contributions including lubrication forces are included. See Refs. 30, 44, and 48 for details on this method.

SANS experiments were conducted on the D-22 and D-33 beamlines at the Institut Laue-Langevin (ILL) in Grenoble, France, as well as the NGB30mSANS instrument at the NIST Center for Neutron Research (NCNR) in Gaithersburg, MD, following previously reported protocols and methods.^4,11 The scattering intensity was obtained for scattering vector magnitudes q ranging from 0.004 Å⁻¹ to about 0.5 Å⁻¹. All samples were held in standard quartz Hellma cells at ILL and custom titanium cells with quartz windows at NCNR. Low concentration samples were studied using cells with a 2 mm path length to enhance intensity, while concentrated samples were studied in cells with a 1 mm path length. All lysozyme concentrations were studied at three temperatures (5 °C, 25 °C, and 50 °C). All raw data files were analyzed using software provided through the NCNR.⁴⁹ The resulting reduced data for lysozyme samples were fitted using an isotropic scattering function.

Lysozyme samples were obtained from MP Biomedicals and subsequently purified to remove impurities and excess counter-ions to minimize the solution ionic strength. Purification was conducted by dialyzing reconstituted lysozyme against deionized water at 4 °C, until the resistance of the water reached approximately 18.0 M ℧. This typically required seven changes of deionized water over the course of 48 h. The purified lysoyzme was then lyophilized by freeze drying. Samples were prepared by dissolving purified lyophilized lysozyme in deuterium oxide (D₂O) at 25 °C and gently vortexing to enhance dissolution and homogenization. Samples were subsequently filtered with 0.22 μm syringe filters to remove additional impurities. Protein content was initially determined by the mass fraction of purified lyophilized lysozyme in deuterium oxide, x_L. The intrinsic volume fraction of lysozyme, ϕ_L, is then calculated according to the specific volume, ν₀, reported in the literature (ν₀ = 0.717 ml/g in Ref. 50), according to

where ρ_D is the mass density of D₂O.

Neutron spin echo experiments were performed on the IN-15 beamline at the ILL in Grenoble, France. Samples were prepared by following the procedure described earlier, then pipetted into 1 mm square quartz cells and stored in a custom temperature controlled sample chamber. All samples were thermally equilibrated for at least 30 min at each of the temperatures studied. For our experiments, the instrument was configured to obtain intermediate scattering functions⁵¹ at correlation times up to 50 ns with 30-35 points for q range between 0.03 Å⁻¹ to 0.20 Å⁻¹ at each sample condition studied.

As previously shown, SALR systems can be in a dispersed fluid state, random percolated state, clustered fluid state, and cluster percolated state.²⁹ We therefore first identify the phase states of the studied seven lysozyme samples. Using the experimental potential parameters obtained previously,^11,31 Monte Carlo computer simulations yield the cluster size distribution, N(s), from which the phase points of the seven samples can be determined. We then investigate here to what extent state-of-the-art theoretical methods developed for spherical colloidal particles are capable of predicting the measured transport properties of lysozyme protein dispersions. In the calculation of the dynamic properties, we model the proteins as spherical with isotropic direct and hydrodynamic interactions, which allows for usage of our semi-analytic calculation methods by which the effects of varied interaction parameters such as the attraction strength can be easily studied. By comparing our theoretical predictions for an isotropic SALR model with experimental measurements on lysozyme solutions, we are in the position, first, to quantify and understand the limitations of the model and, second, to assess the importance of intra- and inter-protein structure (anisotropic shape and surface charge distribution and cluster formation, respectively) of lysozyme with respect to their dynamics. For a clearcut analysis of the theoretical predictions, we intentionally omit the usage of any fit parameter in the calculation of transport properties. The presented results can provide guidance for future theoretical and simulation studies and for future improvements and extensions of the employed theoretical methods.

The equilibrium microstructure of the lysozyme samples addressed in the present paper has been investigated in earlier studies, using small angle neutron scattering (SANS), by part of the present authors.^11,31 The resulting HCDY potential parameters and the volume fractions ϕ, intrinsic volume fractions ϕ_L, and temperature T of the seven considered lysozyme solutions are taken from the previous results and listed in Table I for the convenience of discussions in this paper.

We restrict here our discussion of the S(q) data to features required for the experimental-theoretical comparison of transport properties presented in Sec. III B. Note that S(q) and g(r) constitute the only input to our semi-analytic methods for calculating H(q), D(q), and η.

The experimental structure factors S(q) (symbols) for the lysozyme systems under consideration are plotted in Figs. 1(a) and 1(b) based on the data given by Refs. 11 and 31. Both figures demonstrate the dependence of S(q) on temperature. The experimental S(q)’s of samples 1, 2, and 5 have a distinct IRO peak at q_cσ ≈ 3. An IRO peak is indicative of intermediate range structural correlations of clusters and monomers coexisting in solution, which are hypothesized to influence the short-time dynamics. In contrast, samples 3, 4, and 6 show only a weak IRO peak or shoulder while no clear indication of an IRO peak is observed for sample 7. Although intermediate range order is weak or not present under these conditions, cluster formation is still expected to play a role in the corresponding short-time dynamics.

SANS data of lysozyme samples have been fitted using a two-Yukawa model with the modified reversed hybrid mean spherical approximation (rHMSA) closure to solve the Ornstein-Zernike equation with the implemented thermodynamic self-consistency.^11,31,40 Our structure factor calculations based on the ZH closure are in perfect agreement with those calculated with the modified rHMSA closure. Therefore, we directly take the values of the interaction potential obtained previously as the inputs to calculate the ZH closure structure factors which in turn are the inputs for our dynamical theoretical schemes. Shown in Figs. 1(a) and 1(b) are the theoretical calculations of S(q) (lines) using the parameters in Table I with the ZH scheme combined with the decoupling approximation by assuming spherical particles of hard-core diameter σ ≈ 30.74 Å.⁵² 

For all samples except 2 and 5 (at 5 °C), the agreement between the theoretical curves of S(q) and the experimental data is good at q-values near and below the first peak in S(q), where the calculations are most sensitive to changes in the interaction parameters. However, the calculations shown for samples 2 and 5 still nearly quantitatively capture the observed q-dependence of the experimental data. Larger deviations are visible for all samples only at larger q-values, in the region around the second maximum of S(q) where the statistical error bars in the SANS data are quite large. Additionally, deviations between SANS data and theoretical structure factor curves at intermediate q-values can arise from the structural anisotropy of individual lysozyme proteins at high concentrations for which the decoupling approximation may not be valid any more.

The cluster size distribution (CSD) of the seven samples from our MC simulations are shown in Fig. 2(a). We follow previously developed guidelines to analyze the phase state of a system according to the CSD obtained from MC simulations.²⁹ Samples 1-4 and sample 7 have a monotonically decreasing N(s) that is a characteristic of the monomer-dominated dispersed-fluid phase state where highly transient clusters are present in the system. Samples 5 and 6 are in the random percolated phase state, and their CSD functions show a characteristic peak near the total number of particles, N_p = 1728. Small clusters (10 particles or less) have a fractal dimension of about 1.5, while all larger clusters have a consistent fractal dimension of about 2 (not shown here). Similar fractal dimensions of clusters in SALR systems have been reported.⁵³ 

With regards to the dynamics, it is important to note that the CSD calculated by connectivity does not necessarily represent dynamic clusters that can be measured with neutron spin echo.⁵ For example, simple hard-sphere fluids show a similar monotonically decaying CSD to that found in the dispersed fluid state.²⁹ However, the “clusters” found in hard-sphere systems are merely statistical density fluctuations that dissipate over very short time scales.

The CSD is useful when combined with the nearest neighbor number distribution (NND) function, f(N), which is also included in Fig. 2(b), to identify structural order internal to clusters. Here, f(N) is the average fraction of particles having N neighboring particles whose centers are located inside the attractive well range of the SALR potential. Except for the same at low concentration (1), peaks in f(N) are observed at N > 0 for all samples. In general, as T is lowered and ϕ increased, the NND functions in Fig. 2(b) are shifted to larger values, and under similar conditions, N(s) is shifted to larger cluster sizes. As clusters grow, a larger fraction of particles will be located on the inside of the cluster where neighboring particles are more prevalent. While any non-zero number of neighbors will affect the mobility of individual proteins, the influence will be minimal if the lifetime of those “bonds” is shorter than the diffusive time scale, τ_d.

Both the local and intermediate lengthscale structures, represented by the NND function and IRO peak formation, respectively, will influence the short-time dynamics by affecting the time scale of structural rearrangement. Structural heterogeneity of systems in phase states with strong IRO peaks provides a unique structural feature that can affect the dynamics of a system. One example from prior work is the observation of locally glassy behavior in very concentrated lysozyme samples at low temperature exhibiting IRO, where particle motion is sub-diffusive on the time scale of τ_d but becomes diffusive over long time scales. These dynamic features are rationalized by the structural heterogeneity due to IRO, where locally dense regions are separated by low density voids that allow for largescale rearrangement over long time scales.¹¹ 

Having determined the state diagram locations of the studied samples, the short-time dynamics in the form of H(q) of lysozyme samples are calculated using the S(q)’s discussed in Sec. III A, with focus on the samples in the dispersed fluid state. NSE measures the intermediate scattering function, S(q, t), from which the collective short-time diffusion function, D(q), can be determined from the initial slope of S(q, t).^5,51,54 Then, H(q) can be evaluated using the calculated S(q) by the relationship H(q) = D(q)S(q)/d₀ and compared with theoretical predictions.

The comparison of the NSE-determined H(q) of the lysozyme solutions with our BM-PA theory predictions based on the ZH-calculated g(r) and S(q) as input is presented in Fig. 3. The experimental q-range is constrained here to the IRO peak region. Note that one significant advantage of the NSE method in comparison to dynamic light scattering (DLS) is a larger range of accessible wavenumbers and the ability to probe wavenumbers corresponding to the nearest neighbor and IRO distances.

Our theoretical results for H(q) in Fig. 3 are in semi-quantitative agreement with the experimental data. For sample 1 having the lowest considered volume concentration of ϕ = 0.0398, the agreement is nearly quantitative. Notice in particular the distinct IRO peak in the experimental H(q) at q_cσ ≈ 3 that nicely confirms our earlier theoretical prediction of such a low-q peak of the hydrodynamic function in Ref. 30. Note further that there is no adjustable fitting parameter used in the calculation of H(q). The quantitative agreement indicates that the theoretical BM-PA scheme used here can describe the dynamics of the SALR samples very well at relatively low concentrations.

With increasing ϕ, the theoretical predictions become qualitative, and the deduced H(q) is increasingly overestimated as compared with the experimental results. Yet, the calculated H(q) using BM-PA still reflects the most important trends observed in the experimental data, such as the q-values of the crossing points of the different curves. This holds even for the randomly percolated samples 5 and 6 and the high-temperature sample 7. In recent work by part of the authors³⁶ on a SALR dispersion different from the HCDY system considered here, the good accuracy of the hybrid BM-PA scheme has been shown in comparison with elaborate multi-particle collision dynamics simulations to persist for ϕ values extending at least up to ϕ ≈ 0.1.

In Fig. 4, we compare the NSE results for the experimentally directly obtained short-time diffusion function, D(q), with our theoretical predictions. The agreement between theory and experiment is of similar quality as that for H(q) in Fig. 3. In view of Figs. 1(a), 1(b), and 4, and owing to the small values of S(q) for low q-values (low osmotic compressibilities), the deviations in the H(q)’s between theory and experiment are amplified after the division by S(q).

Note that a source of inaccuracy in inferring d_s from NSE measurements is the limited q-range probed experimentally. However, the oscillations both in the experimental and theoretical D(q)’s are small for q ≳ q_c, with the statistical errors in the NSE-D(q) for larger q being comparable in magnitude to the amplitude of the theoretical D(q) at these wavenumbers. Here, q_c is the wavenumber location of the IRO peak in S(q). In fact, earlier NSE measurements on lysozyme solutions covering larger q-values revealed basically a plateau in D(q) at large wavenumbers.^4,7 Hence, we expect the inaccuracy in inferring d_s from the large-q extrapolation of low-q experimental data using d_s/d₀ = D(q → ∞), to be quite small. This leveling off of D(q) may be facilitated by the averaging over the distribution of the orientations of non-spherical particles (i.e., the orientational polydispersity), which can affect the measured (effective) S(q) and D(q) akin to size polydispersity, causing the damping of large-q oscillations.^8,55

The differences between the theoretical predictions and the experimental results for H(q), visible in Fig. 3, can be attributed to the inherent simplifications in our theoretical model, in which the non-spherical shape and patchiness of the lysozyme proteins have been disregarded. Note that the S(q)’s are overall well-described by the spherical HCDY model even for the larger ϕ values, and even for the two systems considered in the random percolated state, as it is visible in Fig. 1(a) and likewise in Refs. 5, 6, 31, and 35. This is due to the orientational averaging invoked in S(q). In contrast, the protein asphericities have a significantly stronger dynamic influence, which grows with increasing ϕ and decreasing temperature, owing to a stronger coupling of translational and rotational protein motion even for short times. This dynamic translation-rotation coupling slows the collective diffusion on lengthscales 2π/q characterized by H(q), and it significantly contributes to the overestimation of the experimental H(q) by the theoretical predictions, which are based on the isotropic pair potential in Eq. (1) and isotropic hydrodynamic no-slip boundary conditions. In this context, using computer simulations, Bucciarelli et al.⁵⁶ recently revealed a strong slowing effect of directional attractive interactions on the cage-diffusion coefficient, D(q_m) = d₀H(q_m)/S(q_m), with q_m > q_c denoting the principal structure factor peak position in these systems, as compared to purely isotropic attractive interactions of comparable strength. Moreover, Roos et al. in Ref. 57 find experimentally for lysozyme solutions that the reduced rotational self-diffusion coefficient, $dr0/dr$⁠, with $dr0$ denoting the rotational diffusion coefficient, d_r, at infinite dilution, is approximately equal to the inverse of the reduced zero-frequency viscosity, η/η₀, and to the reduced long-time translational self-diffusion coefficient d_L/d₀. This approximate equality is suggestive of strong correlations in the dynamics of neighboring particles due to rotational self-diffusion.⁵⁷ Such a peculiar behavior is not observed in suspensions of charged colloidal particles, where 1/d_r has a weaker ϕ-dependence than η and 1/d_L.⁵⁸ 

The formation of large percolated or non-percolated clusters in the lysozyme solutions studied here at larger ϕ and lower T has an additional (slowing) influence on short-time diffusion properties that is not adequately captured by our simple spherical model where the pair distribution function, g(r), is the only microstructural input. While some additional microstructural information is included in N(s) and f(N) as discussed before, these global equilibrium distribution functions can only provide hints on how diffusion is dynamically influenced by the presence of clusters. An enhanced hindrance of collective diffusion [i.e., an overall smaller H(q)] can be expected to occur with increasing mean cluster size and number of next neighbors. Additionally and quite importantly, however, the shape of H(q) is distinctly influenced by the distribution of cluster shapes and their densities, their intra- and inter-cluster dynamics, and here in particular the distribution of cluster lifetimes. Multiparticle collision (MPC) dynamics simulations in Ref. 36 performed for a SALR dispersion model of spherical Brownian particles interacting by a generalized Lenard-Jones-Yukawa potential show that the smaller clusters change their shape significantly on a time scale comparable to τ_d, indicating that the dynamics of cluster phases cannot generally be described by a polydisperse mixture of rigid clusters. For compact clusters having many next neighbor contacts, genuinely non-pairwise additive many-particle HIs also come into play, which are not adequately described by the BM-PA scheme with the d_s part treated as pairwise additive. A quantitative exploration of the aforementioned dynamic influence of the cluster shapes and dynamics on collective diffusion is a demanding task that should be addressed in a future dynamic simulation study for which the present work can provide guidance.

The zero-frequency (low-shear rate) viscosity is calculated as the sum η = η_∞ + Δη. In contrast to H(q), D(q) and the high-frequency viscosity part η_∞, η is a long-time dynamic property influenced additionally by relaxation (memory) effects embodied here in the shear stress relaxation viscosity contribution Δη. We compare our theoretical predictions for η, obtained using the PA scheme for η_∞ and the simplified mode-coupling theory expression in Eq. (12) for Δη, with earlier measurements of the lysozyme solution viscosity taken from Refs. 11 and 31.

Interestingly, as shown in Fig. 5, our hybrid MCT-PA method results for η are in reasonable agreement with the experimental viscosity measurements, in particular at smaller ϕ. From our PA method theoretical results for η_∞ depicted in the inset of Fig. 5, and in view of Δη = η − η_∞, we see that the increase of η with increasing ϕ is mainly due to the shear relaxation viscosity part Δη. We further note that η_∞ is quite insensitive to changes in the pair potential caused by temperature variations. This is in line with the general observation that Δη is more sensitive to the interaction parameters than η_∞.⁴⁸ As one expects, the viscosity is largest for the lowest considered temperature, and it grows strongly with increasing ϕ in particular for the random percolated cluster systems. This viscosity enhancement is not adequately described by our simple model. Notwithstanding these limitations, the first-order MCT expression in Eq. (12) combined with the PA expression for η_∞ is seen to be a valuable and easy-to-implement tool for the calculation of η, giving results in reasonably good agreement with experimental viscosity data for solutions in the dispersed-fluid phase state at small to intermediately large concentrations.

To explore in more detail the interplay of short-range attraction and long-range repulsion on the structure and dynamics of the considered systems, we focus here on the effect of adding attraction to an initially purely repulsive pair potential. This is accomplished by enlarging the attraction strength, K₁, in the HCDY potential stepwise from zero to its full value given in Table I, while all other parameters including the ranges of repulsion and attraction are kept constant. For samples 1 and 6, the respective theoretical predictions for g(r), S(q) and H(q) are shown in Fig. 6, and the theoretical predictions for η and its high-frequency and shear relaxation parts are included in Fig. 7. As argued earlier, while the theoretical results for the dispersed-fluid sample 1 can be expected to be quantitative, the results for the random percolated system 6 are only qualitative but remain useful for identifying general trends.

Consider first, as shown in the figure parts (a) and (b) of Fig. 6, the changes in the ZH-calculated g(r) induced by increasing stepwise the attraction strength. In the limiting case of zero attraction, in (b) there is a well-developed first-neighbor shell peak visible at r/σ ≈ 0.8 induced by the long-range Yukawa repulsion. Owing to the low concentration, the next-neighbor peak in (a) is quite shallow and situated at a larger pair distance than that in (b). In both cases, there is non-zero relative probability [i.e., a non-zero contact value g(σ⁺) > 0] of finding particle pairs in contact. The contact value of g(r) rises sharply with increasing attraction strength, to the value 14.0 in (a) and 15.0 in (b) when the full strength 1.0 × K₁ is attained. While ϕ ≈ 0.04 in (a) is quite small, the g(r) attained for full attraction strength is distinctly different from the zero-concentration form g₀(r) = exp{−V(r)} (dotted line), which has a contact value of 5.8 in (a) and 10.9 in (b). In the g(r) of the concentrated sample 6 in (b), the first next-neighbor shell maximum at r/σ ≈ 1.3 gradually changes into a minimum with increasing attraction strength, reflecting the buildup of a particle depletion zone next to the attraction-enhanced next-neighbor contact layer. Additionally, a characteristic small peak at r/σ = 2 develops with increasing attraction quantifying the attraction-induced tendency of finding an in-line configuration of three touching particles.

The corresponding changes in S(q) induced by adding short-range attraction are depicted in figure parts (c) and (d), respectively. As shown, an IRO peak of S(q) develops with increasing attraction strength [i.e., a low-q shoulder at qσ ≈ 3 in the case of (b)] that replaces the diminishing principal peak of the purely repulsive reference system with 0 × K₁ that is located at a larger wavenumber. Notice here the attraction-independent isosbestic points of S(q) where the structure factor curves for different attraction values intersect. While in general the undulations in S(q) are reflected in those of H(q) with the latter having smaller amplitudes, owing to the significant decrease of the short-time self-diffusion coefficient d_s with increasing attraction strength [c.f. the large-q values of H(q) in figure parts (e) and (f)], there are no corresponding isosbestic points in the depicted hydrodynamic functions. Moreover, the low-q shoulder of H(q) in (f) at qσ ≈ 3 is less pronounced than the corresponding low-q structure factor shoulder in (d). The formation of (transient) dimers and clusters of several particles with increasing attractions causes the sedimentation coefficient K = H(q → 0) to increase since the hydrodynamic friction by backflowing fluid is reduced. And indeed, a small to moderately pronounced increase of H(0) is observed in (f) and (e), respectively.

According to Figs. 7(a) and 7(b) including our theoretical viscosity predictions for samples 1 and 6, respectively, the PA calculated high-frequency viscosity η_∞ increases only weakly, and monotonically, with increasing attraction strength K₁. In contrast, the MCT-calculated shear-relaxation viscosity part, Δη, decreases mildly at smaller K₁ values, but when K₁ is further enlarged, it becomes significantly enhanced (see insets). The slightly non-monotonic attraction-strength dependence of Δη is reflected in a correspondingly non-monotonic K₁ dependence of the zero-frequency viscosity η (filled circles).

We have assessed the applicability of our semi-analytic methods for calculating dynamic properties and viscosity of globular particle dispersions by comparison to experiments on lysozyme solutions. These methods directly relate bulk diffusion and viscosity to the structure factor without additional inputs. We find a favorable comparison between our theoretical predictions of H(q), D(q), and zero-frequency experimental viscosity data for η with the experimental results. This comparison is made without invoking any fit parameter.

There is a clear quantitative agreement between the BM-PA H(q) and the experimental data observed at the lowest considered concentration. At larger ϕ, qualitative agreement is still maintained even though the quantitative agreement is lost. The experimental-theoretical deviations in D(q) and H(q) can be attributed to the neglect of asphericities in the shape and associated electrostatic interactions of proteins and (transient) protein clusters in the theory, as well as to the disregarded patchiness of the short-range attraction contribution, which for larger concentrations lead to a strong coupling of the rotational and translational particle and cluster dynamics. Our findings are in line with earlier experimental work⁵⁷ on lysozyme solutions, where a d_r ∝ 1/η scaling is observed for the lysozyme solutions, but not for spherical colloidal particles with isotropic interactions. The strong influence of particle surface patchiness on short-time diffusion was highlighted recently in a joint experimental-simulation study of the cage-diffusion coefficient D(q_m) of γ_B-crystallin.⁵⁶ The authors of this study show that patchy attractive interactions give rise to smaller values of D(q_m) than those observed for isotropic attractive interactions of comparable strength, which is consistent with our theoretical prediction for the diffusion being faster than experimentally observed.

Quite unexpectedly, our theoretical predictions for the zero-frequency viscosity η are in reasonably good agreement with rheological data for lysozyme solutions, showing that the employed simplified MCT-PA method is an efficient and easy-to-implement tool for viscosity calculations of globular protein solutions even with pronounced intermediate range order. Lysozyme-specific anisotropic interactions, disregarded in our spherical model, are seen to be more influential on short-time diffusion properties than on the shear viscosity, with the latter including also long-time (shear stress relaxation) contributions.

While being approximate, a virtue of the theoretical model presented here is its analytic simplicity that allows for a straightforward identification of trends in the structure, short-time diffusion, and rheology of globular protein solutions in their dependence on the different interaction and system parameters. We have demonstrated this by considering the effects of a gradual increase of the strength K₁ of the short-range attraction contribution, in particular, regarding H(q) and the solution viscosity η where for the latter a slightly non-monotonic K₁ dependence is predicted.

Finally, a few comments on the expected cluster relaxation dynamics are in order here. In a recent simulation study³⁶ where some of the present authors were involved, the intra-cluster dynamics and its inter-relation with the full solution dynamics are explored to some extent, for SALR-LJY model systems belonging to the cluster-fluid and dispersed-fluid phases. According to this MPC simulation study where the salient HIs are fully accounted for, the reversibly formed clusters undergo pronounced conformational and size changes within a time span of a few τ_textd values. The intra-cluster dynamics and lifetime and the solution dynamics as a whole are significantly influenced by the HIs, which are not screened inside the dynamic clusters. The simulations yield for the cluster-fluid phase systems significantly smaller H(q) values than for the dispersed-fluid phase systems, and a distinctly stronger increase of the mean cluster lifetime with increasing strength of attraction between the particles. While the simulation results in Ref. 36 are for LJY spherical particles only, the same qualitative behavior can be expected for the HCDY systems considered in the present theoretical treatment. Systems in the random percolated phase state, and the effects of non-sphericity and surface patchiness of globular proteins, are not considered in Ref. 36. A clear distinction and comprehensive quantification of particle non-sphericity and patchiness effects and of cluster dynamics effects on protein diffusion and rheology in different equilibrium and non-equilibrium SALR phases have to await future (hydro)-dynamic simulations and theoretical studies where these effects will be individually included and explored. Our plan is to perform such studies in the future and to generate benchmark simulation results for theoretical studies using more detailed non-spherical models.

J.R. acknowledges support by the International Helmholtz Research School on Biophysics and Soft Matter (IHRS BioSoft). G.N. and J.R. thank R.G. Winkler, S. Das (all FZ Jülich), R. Roa (Helmholtz-Zentrum Berlin), and M. Heinen (Universidad de Guanajuato, Leon, Mexico) for many helpful discussions. This manuscript was prepared under the partial support of the cooperative Agreement Nos. 70NANB12H239 and 70NANB10H256 from NIST, U.S. Department of Commerce. Access to NGB30mSANS was provided by the Center for High Resolution Neutron Scattering, a partnership between the National Institute of Standards and Technology and the National Science Foundation under Agreement No. DMR-1508249. Certain commercial equipment, instruments, or materials are identified in this document. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology nor does it imply that the products identified are necessarily the best available for the purpose.