To investigate the cooperativity of hydrophobic interactions, the potential of mean force of two- and three-molecule methane clusters in water was determined by molecular dynamics simulations using two methods: umbrella-sampling with the weighted histogram analysis method and thermodynamic integration. Two water models, TIP3P and TIP4P, were used, while each methane molecule was modeled as a united atom. It was found that the three-body potential of mean force is not additive, i.e., it cannot be calculated as a sum of two-body contributions, but requires an additional three-body cooperative term. The cooperative term, which amounts to only about 10% of the total hydrophobic association free energy, was found to increase the strength of hydrophobic association; this finding differs from the results of earlier Monte Carlo studies with the free energy perturbation method of Rank and Baker (1997). As in the work of Rank and Baker, the solvent contribution to the potential of mean force was found to be well approximated by the molecular surface of two methane molecules. Moreover, we also found that the cooperative term is well represented by the difference between the molecular surface of the three-methane cluster and those of all three pairs of methane molecules. In addition, it was found that, while there is a cooperative contribution to the hydrophobic association free energy albeit a small one, the errors associated with the use of pairwise potentials are comparable to or larger than this contribution.