Computational modeling of diffusion and reaction of chemical species in a three-dimensional (3D) geometry is a fundamental method to understand the mechanisms of synaptic plasticity in dendritic spines. In this protocol, the detailed 3D structure of the dendrites and dendritic spines is modeled with meshes on the software Blender with CellBlender. The synaptic and extrasynaptic regions are defined on the mesh. Next, the synaptic receptor and synaptic anchor molecules are defined with their diffusion constants. Finally, the chemical reactions between synaptic receptors and synaptic anchors are included and the computational model is solved numerically with the software MCell. This method describes the spatiotemporal path of every single molecule in a 3D geometrical structure. Thus, it is very useful to study the trafficking of synaptic receptors in and out of the dendritic spines during the occurrence of synaptic plasticity. A limitation of this method is that the high number of molecules slows the speed of the simulations. Modeling of dendritic spines with this method allows the study of homosynaptic potentiation and depression within single spines and heterosynaptic plasticity between neighbor dendritic spines.