Mode coupling and anharmonicity in a native fluctuating protein are investigated in modal space by projecting the motion along the eigenvectors of the fluctuation correlation matrix. The probability distribution of mode fluctuations is expressed in terms of tensorial Hermite polynomials. Molecular dynamics trajectories of Crambin are generated and used to evaluate the terms of the polynomials and to obtain the modal energies. The energies of a few modes exhibit large deviations from the harmonic energy of kT/2 per mode, resulting from coupling to the surroundings, or to another specific mode or to several other modes. Slowest modes have energies that are below that of the harmonic, and a few fast modes have energies significantly larger than the harmonic. Detailed analysis of the coupling of these modes to others is presented in terms of the lowest order two-mode coupling terms. Finally, the effects of mode coupling on conformational properties of the protein are investigated.