Large scale coarse-grained molecular dynamics simulations of rubber-filler composite systems which are models of rubber compositions have been conducted on the K-computer. Different structures of filler aggregates are modeled by reverse Monte Carlo method from a part of USAXS data of rubber compounds obtained in the large synchrotron radiation facility SPring-8. The difference of experimental stress-strain curves between end-modified SBR and non-modified SBR has been reproduced qualitatively by simulation. Behaviors of large strain region of experimental stressstrain curves can be explained by the difference of maximum elongation of network chains due to the difference of numbers of bonds between fillers in end-modified SBR and non-modified SBR. Difference of small strain region of experimental SS curves between end-modified SBR and non-modified SBR can be reproduced qualitatively by models that filler aggregates are defined. This result suggests that filler aggregates play an important role for behavior of small strain region of stress-strain curves of filled rubbers.