High-fidelity cardiac models using attribute-rich finite element based models have been developed to a very mature stage. However, such finite-element based approaches remain time consuming, which have limited their clinical use. There remains a need for alternative methods for novel cardiac simulation methods capable of high fidelity simulations in clinically relevant time frames. Surrogate models are one approach, which traditionally use a data-driven approach for training, requiring the generation of a sufficiently large number of simulation results as the training dataset. Alternatively, a physics-informed neural network can be trained by minimizing the PDE residuals or energy potentials. However, this approach does not provide a general method to easily using existing finite element models. To address these challenges, we developed a hybrid approach that seamlessly bridged a neural network surrogate model with a differentiable finite element domain representation (NNFE). Given its importance in cardiac simulations, we applied this approach to simulations of the hyperelastic mechanical behavior of ventricular myocardium from recent 3D kinematic constitutive model. We utilized a cuboidal domain and conducted numerical studies of individual myocardium specimens discretized by a finite element mesh and assigned with experimentally obtained myofiber architectures. Both parameterized Dirichlet and Neumann boundary conditions were studied. We developed a second-order Newton optimization method, instead of using a stochastic gradient descent method, to train the neural network efficiently. The resulting trained neural network surrogate model demonstrated excellent agreement with the corresponding “ground truth” finite element solutions over the entire physiological deformation range. More importantly, the NNFE approach provided a significantly decreased computational time for a range of finite element mesh sizes for online predictions. For example, as the finite element mesh size increased from 2744 to 175615 elements, the NNFE computational time increased from 0.1108 s to 0.1393 s, while the “ground truth” FE model increased from 4.541 s to 719.9 s. These results suggest that NNFE run times can be significantly reduced compared with the traditional large-deformation based finite element solution methods. The trade-off is to train the NNFE off-line within a range of anticipated physiological responses. However, training time would only have to be performed once before any number of application uses. Moreover, since the NNFE is an analytical function its computational performance will be amplified when the corresponding problem becomes more complex.