Deformation at tectonic plate boundaries involves coupling between rock deformation, fluid flow and metamorphic reactions, but quantifying this coupling is still elusive. We present a new two-dimensional hydro-mechanical-chemical numerical model and investigate the coupling between heterogeneous rock deformation and metamorphic (de)hydration reactions. Rock deformation consists of linear viscous compressible and power-law viscous shear deformation. Fluid flow follows Darcys law with a Kozeny-Carman type permeability. We consider a closed isothermal system and the reversible (de)hydration reaction: periclase and water yields brucite. In the models, fluid pressure within a circular or elliptical inclusion is initially below the periclase-brucite reaction pressure, and above in the surrounding. Inclusions exhibit a shear viscosity thousand times smaller than for the surrounding, because we assume that periclase-water and brucite regions have different effective viscosities. In models with circular inclusions, solid deformation has a minor impact on the evolution of fluid pressure, porosity and reaction front. Models with elliptical inclusions and far-field shortening generate higher rock pressure inside the inclusion compared to circular inclusions, and show a faster reaction-front propagation. The propagating reaction-front increases the inclusion surface and causes an effective, reaction-induced weakening of the heterogeneous rock. Weakening evolves strongly nonlinear with progressive strain. Distributions of fluid and rock pressure as well as magnitudes and directions of fluid and solid velocities are significantly different. The models mimic basic features of shear zones and plate boundaries and suggest a strong impact of heterogeneous rock deformation on (de)hydration reactions and associated reaction-induced weakening. The applied MATLAB algorithm is provided.