density $\rho$. As $\varepsilon_\mathrm{bend}$ is further increased above $\varepsilon_\mathrm{bend} = 30$, $K_\mathrm{NA}$ decreases. This is because the non-affine relaxation process is constrained due to the large rigidity of $\varepsilon_\mathrm{bend}$. As a result, the total modulus of $K=K_\mathrm{A}-K_\mathrm{NA}$ also presents a non-monotonic behavior, which is demonstrated in Fig. \[fig3\]a. From $\varepsilon_\mathrm{bend}=1$ to $10^2$, $K$ decreases from $K
\simeq 80$ to $60$, which is caused by the reduction of $K_\mathrm{A}$. Moreover, $K$ increases from $K \simeq 60$ to $65$ above $\varepsilon_\mathrm{bend} = 100$, which is caused by the reduction of $K_\mathrm{NA}$. Therefore, the $\varepsilon_\mathrm{bend}$ dependence of the bulk modulus $K$ is determined by the competition between the density reduction and the increase in the bending rigidity.

Further, the shear modulus $G$ is analyzed for $L=50$ and presented Fig. \[fig3\]c. It can be observed that the bending rigidity strongly affects the shear modulus compared to the bulk modulus. Particularly, above $\varepsilon_\text{bend} = 10^2$, both of the affine $G_\mathrm{A}$ and non-affine $G_\mathrm{NA}$ components considerably increase. As the shear deformation is anisotropic and causes deformations of the angles $\theta$ of polymeric chains, its response is expected to be highly affected by the bending rigidity. Interestingly, contrary to the important increases of $G_\mathrm{A}$ and $G_\mathrm{NA}$, the 