It is well-known that proline (Pro) cis-trans isomerization plays a decisive role in the folding and stabilization of proteins. The conformational coupling between isomerization states of different Pro residues in proteins during conformational adaptation processes is not well understood. In the present work, we investigate the coupled cis-trans isomerization of three Pro residues using bradykinin (BK), a partially unstructured nonapeptide hormone, as a model system. We use a recently developed enhanced-sampling molecular dynamics method (ω-bias potential replica exchange molecular dynamics; ωBP-REMD) that allows us to exhaustively sample all combinations of Pro isomer states and obtain converged probability densities of all eight state combinations within 885 ns ωBP-REMD simulations. In agreement with experiment, the all-trans state is seen to be the preferred isomer of zwitterionic aqueous BK. In about a third of its structures, this state presents the characteristic C-terminal β-turn conformation; however, other isomer combinations also contribute significantly to the structural ensemble. Unbiased probabilities can be projected onto the peptide bond dihedral angles of the three Pro residues. This unveils the interdependence of the individual Pro isomerization states, i.e., a possible coupling of the different Pro isomers. The cis/trans equilibrium of a Pro residue can change by up to 2.5 kcal·mol-1, depending on the isomerization state of other Pro residues. For example, for Pro7, the simulations indicate that its cis state becomes favored compared to its trans state when Pro2 is switched from the trans state to the cis state. Our findings demonstrate the efficiency of the ωBP-REMD methodology and suggest that the coupling of Pro isomerization states may play an even more decisive role in larger folded proteins subject to more conformational restraints.