Competing roles of aggregation and interfacial interactions in sustainable protein/cellulose nanocrystal-reinforced soft composites

Abstract

Renewable protein matrix nanocomposites reinforced with high-aspect-ratio cellulose nanocrystals (CNCs) offer promising alternatives to petroleum-based plastics. However, they exhibit mechanical properties far below theoretical predictions, often approaching the Hashin–Shtrikman lower bound despite filler geometries that should approach upper-bound behavior. This discrepancy suggests that microstructural features not captured in standard homogenization approaches dominate the mechanical response. We develop a hierarchical Mori–Tanaka framework that accounts for two competing effects: CNC agglomeration, which diminishes load transfer, and interphase stiffening at CNC–matrix interfaces, which enhances it. Applying this model to soy protein isolate (SPI) composites with unmodified and polydopamine-modified CNCs, we demonstrate attenuation of the high agglomeration inherent to SPI/CNC composites without diminishing favorable interfacial effects. Phase maps reveal conditions that could shift composite performance toward the upper bound, making SPI/CNC bio-nanocomposites a potential sustainable alternative to petroleum-based plastics.

---

1 Introduction

The emergence of bio-based polymer nanocomposites as potential alternatives to petroleum-derived plastics has created new challenges for composite mechanics. Driven by the environmental persistence of synthetic polymers and growing demand for biodegradable materials,^1–4 protein–matrix composites may offer solutions for certain low-strength engineering applications.^5–8 Among promising sustainable options, soy protein isolate (SPI), derived from soybeans, can be polymerized into film-grade bioplastics,⁹ but deleterious effects of the needed plasticizers like glycerol must still be countered with reinforcement.¹⁰

Cellulose nanocrystals (CNCs), produced through acid hydrolysis of plant-based agricultural residues such as soy leaves and stalks, offer a potential sustainable solution, as they are strong, naturally abundant, and biodegradable.^11–14 However, CNC-reinforced biopolymers systematically underperform theoretical predictions by an order of magnitude or more. CNCs possess aspect ratios exceeding 20 and axial elastic moduli of 50–200 GPa^15,16 properties that, in classical composite theory, should yield effective stiffnesses approaching the Hashin–Shtrikman (H–S) upper bound. Instead, these materials cluster near the lower bound, exhibiting behavior characteristic of compliant particulate inclusions rather than stiff, high-aspect-ratio fibers.¹⁷ Three features may contribute to this discrepancy. First, the modulus contrast between filler and matrix is extreme. While conventional glass-fiber/epoxy composites exhibit ratios of order 10² : 1, CNC/biopolymer systems routinely exceed 10⁴ : 1.^18–20 At high contrast, load transfer becomes sensitive to interface quality and filler connectivity, and small deviations from ideal dispersion produce large mechanical penalties. Second, aggregation is not incidental but intrinsic. The same hydroxyl-rich CNC surfaces that enable favorable matrix interactions also drive strong inter-particle attraction, leading to aggregation even at low volume fractions.^21–23 Surface modification can modulate this, but rarely eliminates it.^24–27 Third, the matrices themselves are sensitive to surface interactions: biopolymers such as SPI form heterogeneous networks with spatially varying stiffness,²⁸ and their interactions with CNC surfaces involve electrostatic, van der Waals, and hydrogen-bonding contributions that are difficult to characterize independently.²⁹

Classical micromechanical models such as the Halpin–Tsai, Mori–Tanaka, and self-consistent model are effective for systems where dispersion can be controlled and interfaces are relatively simple.^30–33 These models assume uniform filler distribution within a homogeneous matrix, treating CNCs as isolated, perfectly bonded inclusions. Such assumptions, when applied to a bio-nanocomposite, predict composite stiffness far exceeding measured values. We address this by developing a hierarchical Mori–Tanaka model that accounts for two competing effects: filler agglomeration, which concentrates CNCs into mechanically inefficient clusters, and interphase stiffening, which enhances load transfer at CNC–matrix boundaries. The model partitions CNCs between dispersed and agglomerated states, characterized by parameters ξ (the volume fraction occupied by agglomerated domains) and ζ (the fraction of total CNC volume within those domains). Interphase effects enter through ρ_inter (the interphase volume relative to agglomerate volume) and a stiffening factor q_s that couples interphase modulus to local CNC concentration. This four-parameter framework captures microstructural features that single-step homogenization cannot represent, while remaining analytically tractable and physically interpretable.

We apply this framework to SPI/CNC nanocomposites, using recently published experimental data as both motivation and validation.¹⁷ In that work, surface modification of CNCs with polydopamine (polyDOPA) increased composite stiffness by more than 50% relative to unmodified CNCs at equivalent loading, while transmission electron microscopy revealed improved dispersion in modified systems (Fig. 1). This compatibilizing effect arises from the multifunctional catechol/quinone and amine groups introduced by polyDOPA.²⁵ These groups enhance physicochemical interactions with the SPI matrix through hydrogen bonding and π–π interactions, while also altering the CNC surface characteristics to improve interfacial affinity and nanoscale adhesion, thereby improving CNC dispersion and interfacial load transfer. These observations suggest that polyDOPA enhances properties through both reduced aggregation and strengthened interfaces, but the experiments cannot separate these contributions.

Fig. 1 Mechanical characterization of cellulose nanocrystal (CNC) reinforced soy protein isolate (SPI)–glycerol composites. (a) Chemical structure and transmission electron microscope (TEM) images showing how CNCs change with surface modification: (i) unmodified CNC regions appear as long, needle-like structures, while (ii) after coating with polydopamine (polyDOPA), CNC regions appear shorter and thicker. Dotted lines represent hydrogen bonds between molecules. (b) Stress–strain measurements for SPI/CNC nanocomposites. Each curve represents a different concentration (weight percentage, wt%) of CNCs. Higher CNC content makes the material stiffer but reduces the failure strain. (c) Summary of material stiffness (Young's modulus) for different CNC concentrations, comparing unmodified CNCs (SPI/CNC) to surface-modified CNCs (SPI/PolyDOPA-CNC). Surface modification with polyDOPA significantly increases material stiffness. Figure panels reproduced with permission from ref. 17.

Our model provides this decomposition. By fitting agglomeration and interphase parameters to the two data sets (unmodified and polyDOPA-modified CNCs) we quantify the relative importance of each mechanism. Phase maps constructed over the parameter space identify conditions under which bio-nanocomposites could approach their theoretical stiffness limits, providing specific targets for matrix modification and processing optimization.

Beyond the immediate application to SPI/CNC systems, this work establishes a modeling approach applicable to the broader class of soft biopolymer nanocomposites. As interest in sustainable materials accelerates, predictive tools that connect microstructure to mechanical performance become essential for rational design. Because soy protein isolate and cellulose nanocrystals are renewable, bio-derived materials, SPI/CNC soft composites may serve as sustainable platforms for applications such as filtration membranes and water purification.^34–36 More broadly, the framework developed here may also be extended to other processing formats, including predicting microstructure–property relationships in 3D-printed bio-nanocomposites, offering a path toward physically grounded prediction that supplements empirical fitting.

2 Modeling

The development of predictive models for CNC-reinforced bio-nanocomposites requires addressing three key challenges identified through study of their experimental characterization:¹⁷ (1) significant CNC aggregation even at low weight fractions, (2) complex interactions at CNC–matrix interfaces that are modified by surface treatments, and (3) a substantial modulus mismatch between CNCs and the SPI–glycerol matrix that affects load transfer. To address these challenges, we first developed a model for ideally dispersed CNCs to establish theoretical bounds, then incorporated effects of aggregation, and finally accounted for interphase regions reported in experimental observations of polyDOPA-modified CNCs.

2.1 Stiffness of nanocomposites with fully dispersed CNCs

Although TEM analysis of SPI/CNC composites reveals significant CNC agglomeration even at 5 wt% CNC loading,¹⁷ surface treatment enhances dispersion within the matrix.³⁰ We first considered an idealized scenario where high-aspect-ratio CNCs are uniformly dispersed within the matrix (Fig. 2). As appropriate for intermediate, sub-percolation CNC volume fractions,^37–39 we followed the Mori–Tanaka approach, where the representative volume element (RVE) consists of a continuum with CNC inclusions of volume fraction V_c perfectly bonded to a homogeneous SPI–glycerol matrix of volume fraction V_m = 1 − V_c. In this method, each inclusion experiences the matrix strain as its far-field strain rather than the overall composite strain. The ensemble average of the strain in the matrix phase acts as the reference strain, and the strain inside each inclusion is related to the matrix strain through the dilute strain concentration tensor. The effective properties are then determined by volume averaging the fields over all phases, considering their interactions through the matrix mean-field. The Mori–Tanaka effective elastic stiffness tensor of the nanocomposite, C̄, is:

C̄ = C_m + V_c (C_c − C_m): A_c. (1)

where bold symbols represent second or fourth-order tensors; the operator “:” indicates an inner product; C_m and C_c are the elastic moduli of the matrix and CNCs, respectively; and the concentration tensor A_c relates ε₀ to the average strain in CNCs ε_c, as ε_c = A_cε₀. A_c is expressed as:

A_c = [I − S: C_m⁻¹: (C_m − C_c)]⁻¹. (2)

where I is the identity tensor, and S is the Eshelby tensor,⁴⁰ which depends on the inclusion's shape and orientation. Therefore, C̄ can be estimated from the stiffness tensors of the two phases (C_m,C_c) and S for the specific inclusion geometry and orientation.

Fig. 2 Schematic of an idealized model for CNC reinforcement in a SPI–glycerol matrix. The model for the idealized case assumes perfectly dispersed CNCs (unmodified or polyDOPA-modified) within a continuous matrix phase. This idealized case serves as a theoretical baseline for understanding the effects of CNC surface modification and evaluating the impact of agglomeration and interface quality on composite mechanics. PolyDOPA-modified CNCs are shown with their characteristic surface coating, which affects both interfacial interactions and dispersion behavior.

Although transverse properties of even highly anisotropic long, slender fibers do not typically contribute strongly to the mechanics of composites,^20,41 we explored the possibility that transverse isotropy of CNCs or CNC bundles might influence the elastic modulus; this effect proved negligible for SPI/CNC nanocomposites as well, even when the transverse elastic modulus differed by an order of magnitude (Fig. 8). For this overly general case, with the 1-direction aligned along the axis of the CNC, the elastic stiffness matrix C_c can be expressed using Hill's notation^42–44 as:

where n_c = E¹¹_c (1 − ν₂₃)/(1 − 2ν₁₂ν₂₁ − ν₂₃), k_c = E¹¹_c/[2(1 − 2ν₁₂ν₂₁ − ν₂₃)], l_c = 2ν₁₂k_c, p_c = E¹¹_c/[2(1 + ν₁₂)], and m_c = E²²_c/[2(1 + ν₂₃)]. For unmodified CNCs, E¹¹_c = 50–200 GPa, and E²²_c = E³³_c = 0.1E¹¹_c.¹⁵ We set ν₁₂ = ν₂₃ = ν₁₃ = 0.25, well below the thermodynamic limits.

For CNCs coated with a layer of polyDOPA (10–15 nm thickness), we used the Voigt model to estimate the effective moduli, Eⁱⁱ_polyCNC = Eⁱⁱ_polyV_poly + Eⁱⁱ_c (1 − V_poly). Here, i = 1 or 2. Eⁱⁱ_poly = 2.5 GPa is the isotropic elastic modulus of the polyDOPA layer.⁴⁵ V_poly = 0.8 is the volume fraction of the polyDOPA coating, which was estimated from published TEM images.¹⁷

For randomly dispersed and oriented CNCs, the bulk modulus K and shear modulus G of the nanocomposite were written following the approach of Benveniste,⁴⁶

where K_m and G_m are the bulk modulus and shear modulus of the SPI–glycerol matrix, and α_c, β_c, δ_c, and η_c are,⁴⁷

Finally, the Young's modulus E of the nanocomposite was calculated as,

Although the SPI–glycerol matrix can form biphasic nanostructures with intertwining SPI-rich and glycerol-rich phases ranging in size from 1 to 5 nm, as observed in previous TEM analysis,²⁸ the diameters of the CNCs, which range from 10 to 20 nm for unmodified CNCs and 30 to 50 nm for polyDOPA–CNCs, are significantly larger than the intertwining phase size. Therefore, we treated the SPI–glycerol matrix as isotropic, with Young's modulus E_m = 14.7 MPa and the Poisson ratio ν_m = 0.45.¹⁷

2.2 Mechanics of nanocomposites with agglomerated CNCs

TEM imaging reveals non-uniform dispersion of CNCs within the matrix: certain local regions contain “agglomerated phases” or CNC “inclusions,” while other regions show dispersed CNCs.¹⁷ The CNC fibril bundles are randomly oriented, with an aspect ratio of approximately 2. To capture this, we used a two-step Mori–Tanaka model based on an RVE containing both agglomerated CNC bundles and dispersed CNCs within the SPI–glycerol matrix (Fig. 3). This hierarchical approach first determined the homogenized stiffness of the agglomerated phase and matrix separately, then calculated the homogenized stiffness of the complete nanocomposite using these effective properties. We used the two non-dimensional parameters of Shi, et al.:⁴⁷where W^aggl_c and W^disp_c are the volumes of CNCs in the agglomerated phases and dispersed phases, respectively, W_aggl is the volume of the inclusions, and W is the total volume of the nanocomposite. Based on eqn (4) and (5), the bulk and shear moduli of the agglomerated phase and medium are:

Fig. 3 Two-step micromechanical modeling approach for CNC-reinforced SPI–glycerol composites with agglomeration. (a) Schematic of composite microstructure, showing CNCs partitioned into two states: agglomerated bundles and individually dispersed CNCs within the SPI/glycerol matrix. The interphase region between CNCs and matrix is highlighted, representing the zone of modified matrix properties due to CNC–matrix interactions. (b) The two-step homogenization process: step 1 determined effective properties of the agglomerated CNC phase and the matrix containing dispersed CNCs separately. Step 2 combined these components using Mori–Tanaka theory to predict overall composite properties. This hierarchical approach enabled explicit consideration of both agglomeration effects and interphase contributions to composite stiffness.

Finally, the elastic stiffness of the nanocomposite, considering agglomeration effect, was estimated as:

Here, ψ and ϕ represent scalar quantities of the Wu tensor T.⁴⁸ For spherical inclusions,⁴⁸

Where

2.3 Interphase stiffening between CNCs and the matrix

A transitional interphase region between CNCs and matrix can enhance stress transfer and increase the effective matrix stiffness.⁴⁹ While AFM indentation experiments have revealed gradients in elastic properties within these interphases,^50,51 we adopted two simplifying assumptions: the interphase was treated as homogeneous and isotropic. These assumptions, while introducing some error, enabled analytical treatment of the problem.

For an interphase volume fraction relative to the agglomerated phase of ρ_inter = W_inter/W_aggl, where W_inter is the interphase volume, rule of mixtures (Voigt) estimates the effective bulk modulus K^s_medium and shear modulus G^s_medium of the medium were written:

K^s_medium = K_medium(1 − ξ − ρ_interξ) + K_interρ_interξ, (20)

G^s_medium = G_medium(1 − ξ − ρ_interξ) + G_interρ_interξ. (21)

The moduli of the matrix within the interphase increase with the degree of aggregation,^51–53 which was taken as proportional to the volume fraction of CNCs within the agglomerated phase (W^aggl_c/W or V_cζ). To model this, we introduced a strengthening factor q_s, representing the stiffening of the interphase due to strong CNC aggregation, so that the bulk modulus K_inter and shear modulus G_inter were written:

K_inter = q_sV_cζK_m, (22)

G_inter = q_sV_cζG_m. (23)

By substituting ξ (1 − ρ_inter) for ξ in eqn (14) and (15), and K^s_medium and G^s_medium for K_medium and G_medium in eqn (16) and (17), we estimated the effects of both domain agglomeration and interphase stiffening on nanocomposite mechanics.

3 Results and discussion

3.1 Comparison with homogenization bounds and estimates

We began by evaluating experimental measurements against theoretical bounds for composite stiffness. The Hashin–Shtrikman (H–S) bounds⁵⁴ provide theoretical limits for the elastic properties of composites with dispersed isotropic, homogeneous inclusions (see complete formulation in Appendix B). Experimental data fell near the H–S lower bound (Fig. 4), with polyDOPA-modified CNCs yielding slightly higher stiffness despite the relatively compliant polyDOPA coating.

Fig. 4 Comparison between experimental measurements and theoretical predictions for the elastic modulus of CNC-reinforced SPI–glycerol composites. Results shown for (a) unmodified CNCs and (b) polyDOPA-modified CNCs. The Hashin–Shtrikman (H–S) bounds represent theoretical upper and lower limits for composite stiffness. The Mori–Tanaka (M–T) predictions assume perfect CNC dispersion. Experimental data (points) fall near the H–S lower bound, well below M–T predictions, suggesting the importance of microstructural features like CNC agglomeration and imperfect interfaces.

This proximity to the lower bound is typical for particulate composites,³⁹ but is unexpected for high-aspect ratio fibers with high stiffness contrast between phases (approximately 10 000 : 1). For well-dispersed, high-aspect-ratio inclusions, properties close to the H–S upper bound are expected.^18,19 Indeed, our Mori–Tanaka estimate for ideally dispersed CNCs (red solid line in Fig. 4) predicted stiffness an order of magnitude higher than that reported experimentally.¹⁷ We explored two potential causes: CNC agglomeration and imperfect interfaces. We therefore examined how these features affect composite stiffness.

3.2 Effects of CNC agglomeration

To study how domain agglomeration affects composite stiffness, we applied our model for several values of ζ (Fig. 5). TEM analysis¹⁷ provided estimates of the agglomeration parameter ξ (0.081 for unmodified CNCs and 0.068 for polyDOPA–CNCs). ξ was estimated from binarized TEM images as the area fraction of the aggregated phase, while ζ represents the fraction of CNCs located within the aggregated domains.¹⁷

Fig. 5 Effect of CNC agglomeration on composite stiffness for (a) unmodified CNCs and (b) polyDOPA-modified CNCs. Curves represent Mori–Tanaka (M–T) predictions for different values of the agglomeration parameter ζ, which denotes the ratio of the volume of CNCs in the agglomerated phase to the total volume of CNCs in the composite. Hashin–Shtrikman (H–S) bounds are shown for reference. Higher ζ values indicate increased agglomeration and correspond to lower composite stiffness. Experimental data points are best fitted with ζ = 0.98 for unmodified CNCs and ζ = 0.88 for polyDOPA-modified CNCs, with agglomeration parameters ξ = 0.081 and 0.068 respectively, where ξ represents the ratio of the volume of the agglomerated phase to the volume of the entire composite. The lower ζ value for modified CNCs indicates improved dispersion, demonstrating the beneficial effect of polyDOPA surface treatment.

As expected, increasing ζ reduced composite stiffness. We then asked what values of ζ would be required to fit experimental data for the two CNC types. ζ = 0.98 was required for unmodified CNCs, while ζ = 0.88 was required for modified CNCs. This difference aligned with the expected effect of surface modification, with improved CNC–matrix affinity leading to better dispersion.

However, the required value of ζ = 0.98 for unmodified CNCs exceeds experimental measurements (ζ_exp = 0.909 ± 0.033 at 10 wt% CNCs).¹⁷ This discrepancy suggests our model may be compensating for other effects. Two possible explanations are: (1) CNC degradation during sonication, which can reduce crystallinity and introduce defects,^55,56 and (2) matrix inhomogeneities, particularly variations in protein crystallinity that could locally reduce matrix stiffness. These effects would require artificially high values of ζ in the model to match experimental stiffness measurements.

3.3 Role of interphase properties

To study how interphase properties affect composite stiffness, we examined the effects of two key parameters: the interphase volume ratio ρ_inter and the stiffening factor q_s. Varying ρ_inter (Fig. 6(a_i) and (a_ii)) revealed that larger interphase regions increased composite stiffness for both unmodified and modified CNCs. Similarly, increasing the stiffening factor q_s at fixed ρ_inter = 0.8 (Fig. 6(b_i) and (b_ii)) enhanced composite stiffness, with polyDOPA-modified CNCs showing greater sensitivity to this effect.

Fig. 6 Mori–Tanaka micromechanical modeling predictions incorporating both agglomeration and interphase effects. (a) Impact of interphase volume ratio (ρ_inter) on composite stiffness for (a_i) unmodified CNCs and (a_ii) modified polyDOPA–CNCs, with fixed strengthening factor q_s = 2000. (b) Effect of interphase strengthening factor (q_s) on composite stiffness for (b_i) unmodified CNCs and (b_ii) modified polyDOPA–CNCs, with fixed interphase volume ratio ρ_inter = 0.8. Hashin–Shtrikman (H–S) bounds are included for reference. Agglomeration parameters are fixed at ξ = 0.081, ζ = 0.98 for unmodified CNCs and ξ = 0.068, ζ = 0.88 for polyDOPA–CNCs. Increasing either ρ_inter or q_s enhanced composite stiffness, with polyDOPA-modified CNCs showing greater sensitivity to these interphase parameters.

These results suggested that surface modification works through two mechanisms: enhancing interfacial interactions and increasing the effective interphase volume. The combined effect can shift composite properties significantly closer to the H–S upper bound (Fig. 6(b_ii)). This finding has important implications for bio-nanocomposite design. Strategic modification of CNC surfaces could enable substantial improvements in mechanical performance without requiring changes to bulk material composition. The following section explores specific strategies for implementing these insights in bio-nanocomposites with better mechanical stiffness.

3.4 Guidelines for bio-nanocomposite design

3.4.1 Matrix and CNC stiffness optimization. The performance of bio-nanocomposites depends strongly on the stiffness contrast between the matrix (E_m) and CNCs (E_CNC). To understand this relationship quantitatively, we examined the ratio of predicted composite stiffness to the theoretical H–S upper bound, E^Upper_H–S/E^aggl_M–T (Fig. 7(a_i) and (a_ii)). This ratio provides a metric for how close a composite comes to achieving its theoretical maximum stiffness.

Fig. 7 Design maps for optimizing CNC-reinforced SPI nanocomposites. (a) Colormaps showing the ratio of actual to theoretical maximum stiffness (E^Upper_H–S/E^aggl_M–T) as a function of matrix modulus (E_m) and CNC modulus (E_CNC) for (a_i) unmodified CNCs and (a_ii) modified CNCs. Contour lines indicate modulus ratios of E_CNC/E_m = 10, 100, 1000, and 10 000. Parameters fixed at V_c = 0.1, with ξ = 0.081, ζ = 0.98 for unmodified CNCs and ξ = 0.068, ζ = 0.88 for modified CNCs. (b) Colormaps showing the influence of interphase volume ratio (ρ_inter) and agglomeration factor (ζ) on normalized modulus for (b_i) unmodified CNCs (ξ = 0.081) and (b_ii) modified CNCs (ξ = 0.068), with V_c = 0.1 and q_s = 2000. These maps provide quantitative guidance for optimizing composite properties through control of matrix stiffness, CNC properties, and processing conditions.

Our analysis revealed three key findings. First, CNC reinforcement becomes more effective as matrix stiffness increases, with particularly dramatic improvements beyond E_m = 100 MPa. This effect is amplified in composites with polyDOPA-modified CNCs, suggesting surface modification enhances load transfer. Second, matrix stiffness (E_m = 14.7 MPa) in Wang et al.¹⁷ falls well below optimal values, indicating substantial room for improvement through reduced water content, increased crosslinking, or decreased plasticizer concentration.²⁸

Third, increasing CNC stiffness consistently improves composite performance, particularly in combination with stiffer matrices. This suggests two complementary paths to improvement: (1) optimizing acid hydrolysis conditions to reduce CNC defects, and (2) using high-temperature annealing to enhance CNC molecular alignment and eliminate surface defects.⁵⁷ Based on these findings, we recommend targeting matrix stiffness above 100 MPa and CNC stiffness above 100 GPa to approach theoretical performance limits. These targets could be achieved through combined improvements in both matrix formulation and CNC processing.

3.4.2 Strategies for optimizing CNC dispersion and interfaces. Phase maps examining the combined effects of agglomeration (ζ) and interphase volume (ρ_inter) reveal two key pathways for improving composite performance (Fig. 7(b_i) and (b_ii)). First, increased agglomeration (ζ) consistently reduces load-bearing efficiency, highlighting the critical importance of achieving uniform CNC dispersion. While conventional mechanical stirring and sonication methods can achieve dispersion, they risk introducing CNC defects. Alternative approaches such as high-pressure homogenization⁵⁸ and surfactant-assisted processing⁵⁹ may provide better dispersion while preserving CNC integrity.

Second, increasing the interphase volume ratio (ρ_inter) enhances composite stiffness, particularly when combined with good dispersion (low ζ). This suggests that surface modifications promoting strong CNC–matrix interactions could significantly improve performance. Chemical approaches such as grafting and functionalization²² offer promising routes to optimize these interfacial regions. The synergistic relationship between dispersion and interface quality indicates that processing strategies should target both features simultaneously rather than addressing them independently.

3.5 Model limitations

The SPI/CNC composite system was chosen as a representative sustainable composite with a renewable protein matrix, bio-derived nanofillers, and tunable interfaces. However, the modified Mori–Tanaka framework can be extended to other soft nanocomposites where filler aggregation and interfacial load transfer govern reinforcement.

While our micromechanics framework captures key features of CNC-reinforced composites, several important simplifications merit attention in future work. Matrix heterogeneity remains unaddressed. The current model treats the SPI–glycerol matrix as homogeneous, neglecting variations in protein structure and crystallinity that affect both load transfer and interfacial interactions.⁶⁰ An iterative Mori–Tanaka approach incorporating these features could provide more accurate predictions, particularly at high CNC concentrations where matrix structure may be important.

Additionally, the parameters associated with concentration dependence are difficult to characterize. ξ was estimated from TEM analysis at 10 wt% CNCs, but agglomeration behavior likely varies with concentration. Systematic TEM studies across multiple concentrations could establish this relationship, enabling more accurate modeling across the full range of practical CNC loadings.

The present modified Mori–Tanaka model assumes dispersed or locally aggregated CNC inclusions embedded in a continuous SPI matrix. Therefore, once CNCs reach the percolation threshold and form a connected load-bearing network, the composite behavior would no longer be captured by the current homogenization framework. If high volume fraction SPI/CNC composites prove useful technologically, future extensions could incorporate effects of percolation.

Finally, the interface physics model here is a simplified representation, neglecting both property gradients and potential damage evolution. More sophisticated approaches incorporating multi-layer or gradient-based models could better capture interfacial complexity. Additionally, integrating cohesive zone models^61,62 might enable prediction of progressive damage and failure, particularly important for understanding composite durability under cyclic loading.

4 Conclusions

This study provides insight into the mechanics of CNC-reinforced bio-nanocomposites through integrated experimental characterization and micromechanics modeling. Surface modification of CNCs with polyDOPA enhanced composite stiffness through two mechanisms: improved dispersion and stronger interfacial interactions. However, composite properties remained closer to theoretical lower bounds than expected for high-aspect-ratio reinforcements, pointing to opportunities for further optimization.

Our micromechanics framework revealed roles of both CNC aggregation and interphase properties. The model successfully captured experimental trends and provided quantitative targets for improvement: matrix stiffness above 100 MPa and CNC stiffness exceeding 100 GPa could enable order-of-magnitude improvements in composite performance. The analysis highlighted specific pathways for optimization: enhanced processing methods like high-pressure homogenization to achieve better CNC dispersion while minimizing defects, strategic surface modifications to control interface properties, and matrix modifications to increase stiffness through reduced plasticizer content or increased crosslinking.

These results inform rational design of sustainable bio-nanocomposites. While challenges remain, particularly in simultaneously optimizing multiple competing factors, this work demonstrates that strategic materials processing guided by mechanical modeling can enable substantial improvements in performance. Future work incorporating more sophisticated interface physics and concentration-dependent effects could further refine these design strategies, potentially enabling bio-based composites to rival conventional petroleum-based plastics.

Conflicts of interest

The authors declare no competing interests.

Data availability

The MATLAB codes used to implement the hierarchical Mori–Tanaka micromechanics model in this study are available on GitHub at: https://github.com/jinhanxun/hierarchical-Mori-Tanaka. Additional information can be provided by the corresponding author upon reasonable request.

Appendices

Appendix A: CNC transverse isotropy has negligible effect on composite moduli

Fig. 8 Effect of CNC transverse isotropy on the elastic modulus of CNC-reinforced SPI–glycerol composites. Mori–Tanaka (M–T) predictions are shown for perfect CNC dispersion, considering two transverse isotropy cases: E₂₂ = 0.1E₁₁ and E₂₂ = E₁₁. Results indicate that CNC transverse isotropy has a negligible effect on composite moduli.

Appendix B: Hashin–Shtrikman bounds

The bulk modulus and shear modulus from Hashin–Shtrikman lower bound (LB) and upper bound (UB) estimates are:

Here, K_m and K_c are the bulk moduli of SPI–glycerol matrix and CNCs, respectively; G_m and G_c are the shear moduli of the matrix and CNCs, respectively; and V_m and V_c are the volume ratios of the matrix and CNCs, respectively. Young's modulus was determined as E = 9KG/(3K + G).

Acknowledgements

This work was funded by the NSF through grants OIA-2219142 and DMR 2105150.

Notes and references

1. N. Gowthaman, H. Lim, T. Sreeraj, A. Amalraj and S. Gopi, Advantages of biopolymers over synthetic polymers: social, economic, and environmental aspects, Elsevier, 2021, pp. 351–372.
2. H. Jin, T. Jiao, R. J. Clifton and K.-S. Kim, J. Mech. Phys. Solids, 2022, 164, 104898.
3. T. H. Epps, L. T. Korley, T. Yan, K. L. Beers and T. M. Burt, JACS Au, 2021, 2, 3–11.
4. X. Zhao, Y. Wang, X. Chen, X. Yu, W. Li, S. Zhang, X. Meng, Z.-M. Zhao, T. Dong and A. Anderson, Matter, 2023, 6, 97–127.
5. A. K. Mohanty, F. Wu, R. Mincheva, M. Hakkarainen, J.-M. Raquez, D. F. Mielewski, R. Narayan, A. N. Netravali and M. Misra, Nat. Rev. Methods Primers, 2022, 2, 46.
6. A. Samir, F. H. Ashour, A. A. Hakim and M. Bassyouni, npj Mater. Degrad., 2022, 6, 68.
7. C. Calvino, N. Macke, R. Kato and S. J. Rowan, Prog. Polym. Sci., 2020, 103, 101221.
8. J. Jeon, Z. Wang, H. Li, M. Senanayake, S. V. Pingali, H. Jin, K. Z. Lee, S. V. Subramani, L. Belaygorod and B. Arif, et al., Small, 2026, e06184.
9. F. Song, D.-L. Tang, X.-L. Wang and Y.-Z. Wang, Biomacromolecules, 2011, 12, 3369–3380.
10. D. Trache, M. H. Hussin, M. M. Haafiz and V. K. Thakur, Nanoscale, 2017, 9, 1763–1786.
11. J. Shojaeiarani, D. Bajwa and A. Shirzadifar, Carbohydr. Polym., 2019, 216, 247–259.
12. A. B. Rashid, M. E. Hoque, N. Kabir, F. F. Rifat, H. Ishrak, A. Alqahtani and M. E. Chowdhury, Polymers, 2023, 15, 4070.
13. S. Cao, P. Rathi, X. Wu, D. Ghim, Y.-S. Jun and S. Singamaneni, Adv. Mater., 2021, 33, 2000922.
14. Y. Habibi, L. A. Lucia and O. J. Rojas, Chem. Rev., 2010, 110, 3479–3500.
15. R. J. Moon, A. Martini, J. Nairn, J. Simonsen and J. Youngblood, Chem. Soc. Rev., 2011, 40, 3941–3994.
16. B.-E. Channab, A. E. Idrissi, Y. Essamlali and M. Zahouily, J. Environ. Manage., 2024, 352, 119928.
17. Z. Wang, H. Li, H. Jin, M. Senanayake, S. V. Pingali, W. Goldberg, D. Kobayashi, G. Genin and M. Foston, Polym. Compos., 2025, 46, 16337–16348.
18. G. Weng, Int. J. Eng. Sci., 1984, 22, 845–856.
19. Y. Qiu and G. Weng, Int. J. Eng. Sci., 1990, 28, 1121–1137.
20. G. M. Genin and J. W. Hutchinson, J. Am. Ceram. Soc., 1997, 80, 1245–1255.
21. Y. Chu, Y. Sun, W. Wu and H. Xiao, Carbohydr. Polym., 2020, 250, 116892.
22. C. Gomri, M. Cretin and M. Semsarilar, Carbohydr. Polym., 2022, 294, 119790.
23. Y. Habibi, Chem. Soc. Rev., 2014, 43, 1519–1542.
24. C. Zhang, Y. Ou, W. X. Lei, L. S. Wan, J. Ji and Z. K. Xu, Angew. Chem., Int. Ed., 2016, 55, 3054–3057.
25. J. F. Rocha, L. H. Hasimoto and M. Santhiago, Anal. Bioanal. Chem., 2023, 415, 3799–3816.
26. S. X. Peng, S. Shrestha, Y. Yoo and J. P. Youngblood, Polymer, 2017, 112, 359–368.
27. A. Babaei-Ghazvini, B. Vafakish, R. Patel, K. J. Falua, M. J. Dunlop and B. Acharya, Int. J. Biol. Macromol., 2023, 128834.
28. H. Tian, G. Guo, A. Xiang and W.-H. Zhong, Polym. Test., 2018, 67, 197–204.
29. C. Miao and W. Y. Hamad, Curr. Opin. Solid State Mater. Sci., 2019, 23, 100761.
30. J. R. Capadona, K. Shanmuganathan, D. J. Tyler, S. J. Rowan and C. Weder, Science, 2008, 319, 1370–1374.
31. V. Favier, H. Chanzy and J. Cavaillé, Macromolecules, 1995, 28, 6365–6367.
32. E. C. Demir, A. Benkaddour, D. R. Aldrich, M. T. McDermott, C. I. Kim and C. Ayranci, J. Compos. Mater., 2022, 56, 1591–1604.
33. G. Josefsson, F. Berthold and E. K. Gamstedt, Int. J. Solids Struct., 2014, 51, 945–953.
34. L. Bai, A. Ding, G. Li and H. Liang, Chemosphere, 2022, 308, 136426.
35. M. Gao, Y. Shang, B. Li and H. Du, Green Chem., 2022, 24, 9346–9372.
36. N. A. Gokhale, C. S. R. Vusa and S. Panda, RSC Sustainability, 2024, 2, 1849–1858.
37. T. Mori and K. Tanaka, Acta Metall., 1973, 21, 571–574.
38. F. Saadat, V. Birman, S. Thomopoulos and G. M. Genin, J. Mech. Phys. Solids, 2015, 82, 367–377.
39. G. M. Genin and V. Birman, Int. J. Solids Struct., 2009, 46, 2136–2150.
40. J. D. Eshelby, Proc. R. Soc. London, Ser. A, 1957, 241, 376–396.
41. R. Maurin, P. Davies, N. Baral and C. Baley, Appl. Compos. Mater., 2008, 15, 61–73.
42. R. Hill, J. Mech. Phys. Solids, 1964, 12, 199–212.
43. R. Hill, J. Mech. Phys. Solids, 1964, 12, 213–218.
44. R. Hill, J. Mech. Phys. Solids, 1965, 13, 213–222.
45. H. Li, J. Xi, Y. Zhao and F. Ren, MRS Adv., 2019, 4, 405–412.
46. Y. Benveniste, Mech. Mater., 1987, 6, 147–157.
47. D.-L. Shi, X.-Q. Feng, Y. Y. Huang, K.-C. Hwang and H. Gao, J. Eng. Mater. Technol., 2004, 126, 250–257.
48. T. Wu, Int. J. Solids Struct., 1966, 2, 1–8.
49. A. Sharma, M. Thakur, M. Bhattacharya, T. Mandal and S. Goswami, Biotechnol. Rep., 2019, 21, e00316.
50. A. Pakzad, J. Simonsen and R. S. Yassar, Compos. Sci. Technol., 2012, 72, 314–319.
51. S. K. Rahimi and J. U. Otaigbe, Polymer, 2017, 127, 269–285.
52. C. I. Idumah and C. M. Obele, Surf. Interfaces, 2021, 22, 100879.
53. R. Sinko, X. Qin and S. Keten, MRS Bull., 2015, 40, 340–348.
54. Z. Hashin and S. Shtrikman, J. Mech. Phys. Solids, 1963, 11, 127–140.
55. G. Molnár, D. Rodney, F. Martoïa, P. J. Dumont, Y. Nishiyama, K. Mazeau and L. Orgéas, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 7260–7265.
56. E. C. Emenike, K. O. Iwuozor, O. D. Saliu, J. Ramontja and A. G. Adeniyi, Carbohydr. Polym. Technol. Appl., 2023, 6, 100337.
57. J. F. Matthews, M. Bergenstråhle, G. T. Beckham, M. E. Himmel, M. R. Nimlos, J. W. Brady and M. F. Crowley, J. Phys. Chem. B, 2011, 115, 2155–2166.
58. S.-Y. Lee, S.-J. Chun, I.-A. Kang and J.-Y. Park, J. Ind. Eng. Chem., 2009, 15, 50–55.
59. M. Shalauddin, S. Akhter, W. J. Basirun, N. S. Anuar, O. Akbarzadeh, M. A. Mohammed and M. R. Johan, Measurement, 2022, 194, 110961.
60. X. Sui, T. Zhang and L. Jiang, Annu. Rev. Food Sci. Technol., 2021, 12, 119–147.
61. X. Liu, C. E. Athanasiou, C. López-Pernía, T. Zhu, N. P. Padture, B. W. Sheldon and H. Gao, J. Appl. Mech., 2024, 91, 011003.
62. G. Guo and Y. Zhu, J. Appl. Mech., 2015, 82, 031005.

---