In optical communications, four-dimensional (4D) modulation formats encode information onto the quadrature components of two arbitrary orthogonal states of polarisation of the optical field. These formats have recently regained attention due their potential power efficiency, nonlinearity tolerance, and ultimately to their still unexplored shaping gains. As in the fibre-optic channel the shaping gain is closely related to the nonlinearity tolerance of a given modulation format, predicting the effect of nonlinearity is key to effectively optimise the transmitted constellation. Many analytical models available in the optical communication literature allow, within a first-order perturbation framework, the computation of the average power of the nonlinear interference (NLI) accumulated in coherent fibre-optic transmission systems. However, all current models only operate under the assumption of a transmitted polarisation-multiplexed, two-dimensional (PM-2D) modulation format. PM-2D formats represent a limited subset of the possible dual-polarisation 4D formats, namely, only those where data transmitted on each polarisation channel are mutually independent and identically distributed. This document presents a step-by-step mathematical derivation of the extension of existing NLI models to the class of arbitrary dual-polarisation 4D modulation formats. In particular, the methodology adopted follows the one of the popular enhanced Gaussian noise model, albeit dropping most assumptions on the geometry and statistic of the transmitted 4D modulation format. The resulting expressions show that, whilst in the PM-2D case the NLI power depends only on different statistical high-order moments of each polarisation component, for a general 4D constellation also several others cross-polarisation correlations need to be taken into account.