Finite element (FE) models of rat skull bone samples were developed by reconstructing the three-dimensional geometry of microCT images and voxel-based hexahedral meshes. An optimization-based material identification method was developed to obtain the most favorable material property parameters by minimizing differences in three-point bending test responses between experimental and simulation results. An anisotropic Kriging model and sequential quadratic programming, in conjunction with Latin Hypercube Sampling (LHS), are utilized to minimize the disparity between the experimental and FE model predicted force-deflection curves. A selected number of material parameters, namely Young's modulus, yield stress, tangent modulus, and failure strain, are varied iteratively using the proposed optimization scheme until the assessment index 'F', the objective function comparing simulation and experimental force-deflection curves through least squares, is minimized. Results show that through the application of this method, the optimized models' force-deflection curves are closely in accordance with the measured data. The average differences between the experimental and simulation data are around 0.378 N (which was 3.3% of the force peak value) and 0.227 N (which was 2.7% of the force peak value) for two different test modes, respectively. The proposed optimization methodology is a potentially useful tool to effectively help establish material parameters. This study represents a preliminary effort in the development and validation of FE models for the rat skull, which may ultimately serve to develop a more biofidelic rat head FE model.