Hyperelastic material simulations are commonly performed in commercial FE codes due to availability of sophisticated algorithms facilitating virtual characterization of such materials in FEA easily. However, the solution time required is longer in FEA. Especially when excitation frequencies do not interfere with structural modes, flexible multibody simulation offers a lucrative and computationally inexpensive alternative. However, it is difficult to directly characterize hyperelastic materials in commercial MBS simulation codes, so the reduced solution time comes at the cost of decreased simulation accuracy, especially if the designer is provided with crude stress - strain test data. Hence, the need is to overcome the drawbacks in FEA and multibody codes, as well as to leverage best of both these codes simultaneously. A methodology is presented where non-linear stiffness properties of the hyperelastic materials are expressed as an analytical function in terms of constants of hyperelastic constitutive material models. The constants are first determined from best fit done to the test data within commercial FE codes. Subsequently, the governing strain energy density functions of the respective hyperelastic materials are differentiated analytically. Lastly, the constants are substituted into the derived nonlinear force-displacement analytical function, which can be directly input to the commercial MBS simulation code. Such a force-displacement analytical model is scalable and can also be combined with viscoelastic simulation. A case study is presented to demonstrate simple use of such nonlinear force - displacement characteristics that can be applied to uniaxial lumped springs in a multibody simulation.