Numerical models for heaving buoy wave energy converters are a fundamental tool for device design and optimization, power production assessment and model-based controller design. Ideally, models are required to be easy to implement, simple, accurate and computationally efficient. Unfortunately, such features are often conflicting and a compromise has to be reached to define an appropriate model structure. A very common choice is to assume a small amplitude of motion and linearize the model. Despite the attractiveness of computational convenience, linear models quickly become inaccurate when large motion occurs. In particular, the implementation of a control strategy, which aims to increase power absorption, enlarges the operational space of the device and significantly enhances the impact of nonlinearities on the model. There are different possibilities to approach the representation of nonlinearities in heaving point absorbers, each of them characterized by a different level of complexity, computational time requirements and accuracy. This paper compares six different methods: one of them fully-nonlinear (implemented in a computational fluid dynamics environment) and the others based on a linear model with the progressive inclusion of nonlinear restoring force, nonlinear Froude-Krylov force and viscous drag.