QSPR Modeling for Mechanical Properties of Reinforced Thermoplastic Starch ()
1. Introduction
Computational investigation applied to polymer chemistry has been on the rise for the last few decades. The computational approach to research, especially in the field of polymer chemistry provides a multitude of benefits compared to the traditional experimental approach. One of the most advantageous attributes provided by computer-aided investigation of polymers is its incredible predictive power [1]. In fact, through various computational techniques, researchers have been able to predict the variation of physiochemical properties of polymers in correlation with the change in their chemical structures: this prediction method is referred to as Quantitative Structure-Property Relationship modelling (QSPR). Given the current environmental context, there is a paradigm shift in the field of materials research; the trend is moving towards biodegradable materials without harmful effects on the environment. Thermoplastic starch is widely regarded as one of the most promising bioplastics for practical applications [2]. Its appeal stems from its low cost, its origin from abundant renewable agricultural resources, and its inherent biodegradability, which collectively position it as a viable alternative to petroleum-based polymers. To date, limited research has been conducted on the quantitative structure–property relationship (QSPR) of TPS-based materials. Although several studies have examined the mechanical behavior of TPS/NR blends experimentally [3]-[7], no prior work has applied quantitative structure–property relationship (QSPR) modeling to predict the mechanical properties of this biopolymer system. This study provides the first demonstration of a monomer-level QSPR framework for TPS-based materials using connectivity-based descriptors generated in the Synthia module. By combining atomistic geometry optimization with descriptor-driven mechanical property prediction, the present work establishes a computational workflow capable of rapidly estimating stiffness, elasticity, and compressibility of TPS/NR blends without the need for iterative laboratory testing. This approach not only fills a significant gap in the computational prediction of biodegradable polymers but also offers a scalable and cost-efficient method that can be extended to other starch-based or rubber-modified biopolymer systems. The methodology proposed here contributes to advancing the rational design of sustainable materials by providing a predictive model that links molecular structure to macroscopic mechanical performance.
2. Computational Theory and Material
2.1. Geometry Optimization Theory
This investigation was performed using Material Studio 2017 Software. This software provides the user with the possibility to manually sketch the molecular structures targeted for analysis using atomistic windows. By using the “forcite” module, the geometry of the newly sketched structures can subsequently be optimized in order to minimize its molecular spotential energy. In fact, there is a significant link between the geometry of a structure and its energy; a stabilized conformation will tend to display a geometry specific to minimal potential energy [8]. Geometry optimization of the target molecules is a recurring process in molecular simulative studies. Optimization with respect to potential energy is always recommended after the sketching of the structure. This is explained by the fact that manual sketching tends to generate excessive energy configurations, which can lead to incorrect calculations. The structural geometry optimization process can be divided into two steps [9]:
Energy evaluation.
Conformation adjustment.
Step 1. Energy evaluation
In this step, the determination and assessment of energy expression is executed for a defined conformation. Energy expressions are generated from structure coordinates in association with a forcefield. This energy expression defines the potential energy surface corresponding to a particular molecular configuration, depending on its atomic coordinates.
To illustrate this notion, let us consider the energy expression for the water molecule; this equation expresses the potential energy surface of a molecular structure of water. The relationship is given by:
[10] (1)
in this example of energy expression, the forcefield when applied, will determine:
The bond lengths (b) and angles (θ).
The functional form (a simple quadratic in both types of coordinates).
The force constants (K).
The reference O-H bond length (b0) and H-O-H angle (θ) are the values for an ideal O-H bond and H-O-H angle at zero energy, which is not necessarily the same as their equilibrium values in a real water molecule.
Step 2: Conformation adjustment
During this phase a progressive conformation rearrangement is achieved to lower the energy expression value. One adjustment may be enough to reach a minimum; otherwise, many thousands of iterations may be required. The required iterations for conformational refinement vary according to factors such as the algorithm’s design, the energy expression employed, and the size of the molecular structure.
2.2. Polymers QSPR Modeling
The Quantitative structure property relationship modeling (QSPR) of polymers makes practical and effective use of mathematical models for rapid property prediction of polymers assumed from their molecular structure. Since QSPR modelling is often implemented on monomer or oligomer level, it enables the quantification of polymers characteristics by using molecular descriptors (such as connectivity indices, topological indices, and other molecular structural features) [11]. Molecular descriptors are computed from the internal composition and arrangement of the polymer subatomic particles. The Synthia module allows quick QSPR modeling of polymer characteristics using systematic observation and experimentation of data. A broad spectrum of properties (thermodynamic, mechanical, and transport properties) can be predicted for large shapeless homopolymers as well as random copolymers [12]. Synthia calculations are executed based on connectivity indices method.
A connectivity index is label that describes a molecule in terms of topology by quantifying how atoms are connected in the molecule. Connectivity indices aim at representing molecular characteristics such as the degree of branching, molecular size or shape, etc. It is an important parameter for the correlative investigation of properties and chemical structure relationship [13].
The theoretical graph-based approach to molecular property analysis begins by constructing the hydrogen-suppressed molecular graph. When dealing with polymer repeat units, additional considerations are required to ensure consistent representation of chain continuation and to prevent truncation errors. To illustrate the procedure, the repeat unit of poly (vinyl fluoride) (PVF) is used as an example. Figure 1 presents this repeat unit along with the corresponding hydrogen-suppressed graph used for calculations.
Figure 1. The repeat unit (a) and corresponding hydrogen-suppressed graph; (b) Indices for the PVF repeat unit.
Next, two atomic indices, δ and δv, are introduced to characterize the bonding and electronic environment of each non-hydrogen atom. The first index, δ, known as the simple connectivity index, corresponds to the number of other non-hydrogen atoms bonded to a particular atom. The second index, δv, incorporates electronic configuration details and is defined by:
(2)
where Zv is the number of valence electrons of the atom, NH is the number of hydrogens attached to it, and Z is its atomic number.
Using these atomic indices, two bond indices, β and βv, can also be defined. They are given by:
. (3)
and
. (4)
Figure 2 summarizes the calculated atomic and bond indices for the PVF repeat unit. In Figure 2(a), the simple connectivity indices δ and β are shown, whereas Figure 2(b) displays the valence-based indices δv and βv.
Figure 2. The atomic and bond connectivity (a); (b) valence indices for the PVF repeat unit.
2.3. Structure-Property Correlations Fomechanical Properties
Correlative calculations executed by Synia module are based on structure-property equations developed by Seitz [14] for mechanical properties of isotropic amorphous polymers (glass transition temperature > 298 K). If only the room-temperature property values are required, equations 2 and 3 can be applied for the approximation of Poisson’s ratio v (298 K) and Bulk modulus B (298 K), respectively. The elasticity modulus (also called Young’s modulus) expresses the material stiffness, or how easily it can be stretched. Low elastic modulus is equivalent to low resistance to stretching or bending, hence high elasticity. The bulk modulus provides insight on the material resistance to the compression force [15]. The higher the numerical value of the modulus, the harder it is to compress the material, and vice versa [16]. Shear modulus is an expression of material rigidity [17]. It defines the material resistance to transverse deformation. It can be used sometimes to gain perspective about elasticity as well [18].
. (5)
. (6)
where:
V (0 K) = Molar volume at absolute zero (in cm3/mole).
Ecohl = Fedor-type cohesive energy (in J/mole).
lm = Length of repeat unit in its fully extended conformation (in cm).
V (T) = molar volume at temperature T (in cm3/mole).
Vw = van der Waals volume of the repeat ut (in cm3).
The results from Equations (5) and (6) can th be inserted into Equation (4) to determine the elastic modulus (E) and shear modulus (G).
(7)
3. Experimental
3.1. Geometry Optimization
Two three-dimensional molecular structures—the TPS repeat unit and the natural rubber monomer—were first drawn manually and then imported into the Forcite module for geometry optimization. Energy minimization was subsequently performed on these structures using 500 optimization cycles, applying the COMPASS force field [19] together with the SMART algorithm. This procedure produced two optimized molecular models, illustrated in Figure 3(a) and Figure 3(b), with hydrogen, oxygen, and carbon atoms represented by white, red, and grey spheres, respectively.
(a) (b)
Figure 3. (a) Optimized molecular model of TPS repeat unit [20]; (b) Optimized molecular model of NR repeat unit (isoprene).
3.2. Mechanical Properties Calculations
Figure 4. Synthia calculation tab layout tab in Material Studio 2017.
For QSPR prediction of TPS/NR blend mechanical properties, optimized structures of TPS and NR were used as input for investigation. Monomer-based QSPR provides accurate trend prediction for amorphous polymer systems when long-range morphology is not the dominant factor in determining mechanical response. The polymer composite made of Thermoplastic starch and Natural rubber was defined from its individual monomeric repeat units. The bulk temperature was maintained at 298 K and molecular weight for the computational bulk was set at 100,000 amu. Subsequently, QSPR predictive calculations of the selected mechanical characteristics (Bulk modulus, Young’s modulus, Shear modulus, and Poisson’s ratio) of the copolymer were performed using “Synthia” module for a broad concentrations range of natural rubber (0.0 to 1.0 mole fraction). The output was a collection of 11 blend systems with various mechanical properties. A basic layout of the synhia calculation tab is illustrated in Figure 4.
4. Results and Discussion
4.1. Energy Minimization Outcome
As shown in the structural geometry optimization reports, the geometry optimization results reveal a significant stabilization of both the thermoplastic starch (TPS) model compound and the natural rubber (NR) structure following energy minimization. For each system, the initial configuration exhibited high internal stresses, as indicated by elevated total energies and large force values. The optimization process successfully adjusted the atomic coordinates to reach low-energy, mechanically relaxed conformations.
1) Thermoplastic Starch (TPS)
The initial TPS structure displayed an extremely high total energy (264768.53 kcal/mol), predominantly arising from non-bonded interactions, especially van der Waals repulsion, which accounted for more than 99% of the initial energy. The very large RMS and maximum force values indicate a highly strained and non-equilibrated starting geometry. After optimization, the total energy dramatically decreased to 19.99 kcal/mol, reflecting the successful relaxation of the structure. Valence contributions (bond, angle, torsion) also reorganized into physically meaningful ranges, with torsional relaxation playing a major stabilizing role (torsion energy dropping from +14 to −21 kcal/mol). The non-bonded energy became modest and largely dominated by electrostatic contributions. The final RMS and maximum force values confirm the attainment of a well-converged minimum. Overall, TPS reaches a stable configuration with a balanced geometry and significantly reduced internal repulsion.
2) Natural Rubber (NR)
The NR chain started from a moderately high initial energy (25.66 kcal/mol), with most of the contribution originating from valence terms—mainly bond stretching and angular distortions. Non-bonded interactions were comparatively low and slightly stabilizing. Following optimization, the total energy decreased to −5.82 kcal/mol, indicating that the minimized NR conformation is thermodynamically favorable relative to the initial structure. The negative torsion term in the final state highlights the role of rotational adjustments around carbon–carbon bonds in reaching an energetically preferred conformation. Electrostatic interactions remained the main stabilizing factor among non-bonded terms. The low RMS and maximum force values demonstrate full convergence and structural relaxation. Both structures underwent profound geometrical relaxation, transitioning from strained, high-energy configurations to stable low-energy conformations.
TPS showed a drastic drop in energy due to the elimination of severe van der Waals repulsions.
NR achieved moderate stabilization, mostly through torsional relaxation and improved valence geometry.
Table 1. Geometry optimization report for TPS and NR repeat units.
Parameter |
TPS (Initial) |
TPS (Final) |
NR (Initial) |
NR (Final) |
Total Energy (kcal/mol) |
264768.53 |
19.99994 |
25.66035 |
−5.82312 |
Valence Energy (kcal/mol) |
1739.37 |
−13.50 |
30.34 |
−0.54 |
Bond |
1602.20 |
0.80 |
27.08 |
0.10 |
Angle |
151.42 |
6.78 |
2.43 |
1.99 |
Torsion |
−14.25 |
−21.07 |
0.83 |
−2.64 |
Cross Terms (kcal/mol) |
−13.70 |
−5.59 |
−1.39 |
−1.57 |
Nonbond Energy (kcal/mol) |
263042.86 |
39.09 |
−3.29 |
−3.71 |
van der Waals |
262917.76 |
4.23 |
2.22 |
1.70 |
Electrostatic |
125.11 |
34.86 |
−5.51 |
−5.41 |
RMS Force (kcal/mol·Å) |
4.22 × 105 |
9.08 × 10−2 |
5.12 × 101 |
9.83 × 10−2 |
Max Force (kcal/mol·Å) |
2.13 × 106 |
2.76 × 10−1 |
1.38 × 102 |
2.53 × 10−1 |
Convergence Status |
Extremely strained; nonbond repulsion dominant |
Fully converged; stable minimum |
Moderately strained; Valence distortions |
Fully converged; stable minimum |
The optimization results (Table 1) show that TPS initially contained severe steric congestion, reflected in extremely high van der Waals repulsion. Optimization reduced the total energy by more than five orders of magnitude, confirming relaxation into a physically meaningful conformation. NR stabilized mainly through torsional relaxation along the carbon backbone. The dramatic reduction in forces and total energy for both polymers confirms successful convergence suitable for subsequent QSPR analysis.
4.2. Mechanical Properties Predictions of TPS/NR Systems
Figure 5. (a) Plot of Bulk modulus variation with Natural rubber content; (b) Plot of Bulk modulus variation with natural rubber content; (c) Plot of Young’s modulus variation with natural rubber content; (d) Plot of Poisson’s ratio variation with natural rubber content.
According to the QSPR correlations generated by Synthia, an overall reduction in rigidity and resistance to compression is observed. This trend is reflected by the decline in both the shear modulus and the bulk modulus as the mole fraction of NR increases (Figure 5(a) and Figure 5(b)).
Furthermore, the elasticity was enhanced by the incorporation of NR to the computational system. This is illustrated by the decrease of the elastic modulus over the raising mole fraction of rubber (Figure 5(c)). The blend bulk became much less brittle because of incorporation of NR to the bulk. The poisson’s ratio variated more towards 0.5 (Figure 5(d)), which is the ideal value for elastic material and soft materials.
The mechanical trends observed in the TPS/NR systems can be rationalized by examining the molecular interactions within the blend. Natural rubber possesses a highly flexible polyisoprene backbone [21], and its incorporation into the TPS matrix introduces segments with significantly higher chain mobility [22]. This flexibility reduces the overall stiffness of the blend, leading to decreases in both Young’s modulus and shear modulus. In addition, NR molecules interfere with the dense hydrogen-bonding network that normally forms between starch chains; this disruption diminishes the cohesive energy density of the TPS phase and reduces its semi-crystalline order. Cohesive energy density (CED) refers to the amount of energy required to separate molecules from each other and is therefore an indicator of intermolecular attraction within a material [23]. In TPS, a high CED reflects strong starch–starch interactions that contribute to stiffness [24]. When NR molecules are incorporated, they disrupt these interactions and lower the overall CED, which in turn reduces the rigidity of the blend. As the hydrogen-bonding density decreases, the polymer matrix becomes more amorphous and exhibits greater free volume, which further lowers resistance to both compression and shear. This effect is reflected in the declining bulk modulus as NR content increases. Simultaneously, Poisson’s ratio approaches values near 0.5—characteristic of rubber-like, nearly incompressible materials—indicating a transition from a comparatively rigid biopolymer to a soft, elastomer-like behavior consistent with classical polymer physics. Together, these molecular-level effects explain the enhanced elasticity and reduced rigidity of TPS/NR blends predicted by the QSPR model.
5. Conclusion
In this work, QSPR modeling of the TPS/NR blend was carried out using individual monomers to represent the full polymer chain. This strategy offers the potential for a rapid and computationally efficient means of predicting the mechanical behavior of TPS-based materials. The model was evaluated by comparing its predictions with existing experimental data. Four mechanical parameters were considered: bulk modulus (compression resistance), Young’s modulus (elastic stiffness), shear modulus (resistance to shearing), and Poisson’s ratio (transverse deformation relative to applied force). The QSPR correlations indicated that incorporating natural rubber into thermoplastic starch enhances elasticity, as shown by the reduction in Young’s modulus with increasing rubber content. At the same time, a decline in rigidity and compression resistance was observed, reflected by decreasing shear and bulk modulus values as the natural-rubber mass fraction increased. These findings are consistent with the experimental results reported by Carvalho et al. (2003) [25], who showed that incorporating NR into TPS increases flexibility and ductility while generally reducing tensile strength and Young’s modulus unless interfacial compatibility is well managed. Their work reported a clear qualitative decrease in Young’s modulus with increasing NR content, and the QSPR predictions in the present study reproduce this same directional trend, further reinforcing the agreement between the model and experimental evidence. Based on these mechanical characteristics, TPS/NR blends have potential applications in biodegradable and single-use products, controlled nutrient-release fertilizers, as well as in adhesives and sealants [26]-[28]. This study provides the first QSPR-based mechanical prediction of TPS/NR biopolymer blends. The approach offers a rapid computational alternative to traditional mechanical testing and may guide material design for biodegradable packaging and controlled-release applications. Future work should incorporate molecular dynamics (MD) simulations to complement the QSPR predictions, allowing deeper insight into chain mobility, hydrogen-bond interactions, and phase behavior at different blend compositions. Additional computational investigations involving compatibilizers or plasticizers could further improve the predictive capability for real industrial formulations. Finally, expanding the dataset to develop machine-learning QSPR models may also enhance prediction accuracy for a wider range of TPS-based biopolymer blends. Despite the advantages of QSPR modeling, several limitations must be acknowledged; QSPR models rely on simplified molecular descriptors, typically based on monomers [29], which prevents them from capturing long-range intermolecular interactions, chain entanglement, and morphological effects such as phase separation [30]. While monomer-based QSPR provides a valuable and efficient screening framework, the final design of materials for specific applications still requires experimental validation or integration with more advanced multi-scale modeling approaches. Moreover, the predictive accuracy of QSPR relationships depends heavily on the quality and relevance of the underlying correlations, and they do not account for temperature-dependent or time-dependent mechanical behavior. These limitations suggest that QSPR predictions should be complemented by molecular dynamics simulations or experimental validation for a more comprehensive understanding of polymer performance.
Acknowledgements
My heartfelt thanks go to my wife, Rosette Baswagha, for her love, patience, and endless support during this journey. Her encouragement has been a constant source of strength. I am equally grateful to my colleagues for their teamwork, valuable input, and helpful discussions, which have played a key role in the success of this work.