ADP-glucose pyrophosphorylase, a key allosteric enzyme involved in higher plant starch biosynthesis, is composed of pairs of large (LS) and small subunits (SS). Current evidence indicates that the two subunit types play distinct roles in enzyme function. The LS is involved in mainly allosteric regulation through its interaction with the catalytic SS. Recently the crystal structure of the SS homotetramer has been solved, but no crystal structure of the native heterotetrameric enzyme is currently available. In this study, we first modeled the three-dimensional structure of the LS to construct the heterotetrameric enzyme. Because the enzyme has a 2-fold symmetry, six different dimeric (either up-down or side-by-side) interactions were possible. Molecular dynamics simulations were carried out for each of these possible dimers. Trajectories obtained from molecular dynamics simulations of each dimer were then analyzed by the molecular mechanics/Poisson-Boltzmann surface area method to identify the most favorable dimers, one for up-down and the other for side-by-side. Computational results combined with site directed mutagenesis and yeast two hybrid experiments suggested that the most favorable heterotetramer is formed by LS-SS (side-by-side), and LS-SS (up-down). We further determined the order of assembly during the heterotetrameric structure formation. First, side-by-side LS-SS dimers form followed by the up-down tetramerization based on the relative binding free energies.