Transit Timing Variations, or TTVs, can be a very efficient way of constraining masses and eccentricities of multi-planet systems. Recent measurements of the TTVs of TRAPPIST-1 led to an estimate of the masses of the planets, enabling an estimate of their densities. A recent TTV analysis using data obtained in the past two years yields a 34% and 13% increase in mass for TRAPPIST-1b and c, respectively. In most studies to date, a Newtonian N-body model is used to fit the masses of the planets, while sometimes general relativity is accounted for. Using the Posidonius N-body code, we show that in the case of the TRAPPIST-1 system, non-Newtonian effects might be also relevant to correctly model the dynamics of the system and the resulting TTVs. In particular, using standard values of the tidal Love number $k_2$ (accounting for the tidal deformation) and the fluid Love number $k_{2f}$ (accounting for the rotational flattening) leads to differences in the TTVs of TRAPPIST-1b and c similar to the differences caused by general relativity. We also show that relaxing the values of tidal Love number $k_2$ and the fluid Love number $k_{2f}$ can lead to TTVs which differ by as much as a few 10~s on a $3-4$-year timescale, which is a potentially observable level. The high values of the Love numbers needed to reach observable levels for the TTVs could be achieved for planets with a liquid ocean, which, if detected, might then be interpreted as a sign that TRAPPIST-1b and TRAPPIST-1c could have a liquid magma ocean. For TRAPPIST-1 and similar systems, the models to fit the TTVs should potentially account for general relativity, for the tidal deformation of the planets, for the rotational deformation of the planets and, to a lesser extent, for the rotational deformation of the star, which would add up to 7x2+1 = 15 additional free parameters in the case of TRAPPIST-1.