We know that two mass particles attract each other with a force
$$F~=~\frac{G M_1 M_2}{r^2}.$$
But what is the reason behind that? Why does this happen?
Gravity – Why is Gravitational Force Proportional to Masses?
gravitymassnewtonian-gravity
Related Solutions
Worse than electrodynamics, general relativity is non-linear, in the sense that the field from multiple sources is not just the sum of fields from each isolated source. Even the simple case you are asking about, which is the two-body problem in general relativity, has not been solved exactly.
An even simpler case is the limit $m_2 \to 0$. In that case only $m_1$ affects the geometry of spacetime, and $m_2$ follows a geodesic in that spacetime. This has been solved exactly. The linked article gives details.
[Addendum] To directly answer the original question for the limit $m_2 \to 0$:
Of course if $m_2=0$, then the force between the two masses is $0$. But the thing about classical gravity is that the acceleration due to gravity is independent of mass. (This is the equivalence principle and is actually one of the starting points of the theory of general relativity.) So it still makes sense to ask what would be the acceleration of a negligible mass (aka, a test body) due to the gravity of another mass. The classical answer is $a_2 = \frac{Gm_1}{r^2}$.
In general relativity, the acceleration of a test body due to the gravity of a single spherical, homogenous, non-rotating mass is given exactly by the Schwarzschild solution, which link you can consult for details. The result of which is that
$$\ddot{r} = -\frac{Gm_1}{r^2} + r\dot{\theta}^2 - \frac{3Gm_1}{c^2}\dot{\theta}^2$$ $$\ddot{\theta} = -\frac{2}{r}\dot{r}\dot{\theta}$$ [CORRECTED AND SIMPLIFIED Jan 2]
where $r$ and $\theta$ are polar co-ordinates centred on the gravitating mass, the dots represent differentiation by the proper time of the test body.
So the first term is just the classical radial acceleration $-\frac{Gm_1}{r^2}$. The terms $r\dot{\theta}^2$ and $-\frac{2}{r}\dot{r}\dot{\theta}$ are the classical centrifugal and Coriolis acceleration for polar co-ordinates.
What's not classical is the extra term $\frac{3Gm_1}{c^2}\dot{\theta}^2$. Finally, there is the fact that differentiation is with respect to proper time of the test body. Different test bodies will experience time differently. They can be related by:
$$(1 - \frac{r_s}{r})\dot{t}^2 - \frac{\dot{r}^2}{(1 - \frac{r_s}{r}) c^2} - \frac{r^2\dot{\theta}^2}{c^2} = 1$$
The constant $r_s = \frac{2Gm_1}{c^2}$ is introduced for simplicity. It is called the Schwarzschild radius of the gravitating mass.
Here the co-ordinate $t$ is introduced as a reference time, so $\dot{t}$ is the rate of change of reference time with respect to proper time of the test body. For a distant ($r \to \infty$), stationary ($\dot{r}=0, \dot{\theta}=0$) test body, this becomes $\dot{t} = 1$, so the reference time can be interpreted as the time measured on a distant, stationary clock.
In the classical case of course every body experiences the same time, but it can also be compared to the special relativistic case, where the equation would be:
$$\dot{t}^2 - \frac{\dot{r}^2}{c^2} - \frac{r^2\dot{\theta}^2}{c^2} = 1$$
So what's new in general relativity is the factor $(1 - \frac{r_s}{r})$. To get some idea of scale, for the Earth, $\frac{r_s}{r}$ is about one and a half parts per billion on the surface of the Earth. (Note the Schwarzschild solution is only applicable outside the gravitating body.)
This is exact only for $m_2 \to 0$, but it remains a very good approximation as long as $m_2$ is much smaller than $m_1$, such as for a planet orbiting a star.
Regardless of the central force type (gravity, spring, contact, electromagnetic) it causes an acceleration of each body proportional to $\frac{1}{m_i}$. As a result the relative acceleration is proportional to $\frac{1}{m_1}+\frac{1}{m_2}$.
When this relative acceleration is used to calculate the effective mass this force sees (to produce the said acceleration) the result is
$$ \boxed{ \mu = \left(\frac{1}{m_1}+\frac{1}{m_2}\right)^{-1} }$$
It is a kind of grouping for masses, the same way resistors are grouped with $\frac{1}{R} = \frac{1}{R_1}+\frac{1}{R_2}+\ldots$, or springs are grouped $\frac{1}{k}=\frac{1}{k_1}+\frac{1}{k_2}+\ldots$. It is all a consequence of the fact they share something common. For us it is a central force, for resistors it is voltage, and for springs it is also force.
When displacements are shared (the same with springs) then the effective mass would be $\mu=m_1+m_2+\ldots$, and if currents are shared the effective resistance would be $R=R_1+R_2+\ldots$.
I hope I have shed some light into your question.
Best Answer
One could explain "well, gravity is the curvature of spacetime due to the mass-energy". But that would only lead to "well, why does mass-energy curve spacetime?" And, should someone produce a proposed answer to that, the follow-up question would have to be "but why is that so?" etc.
At some point though, one must accept that there are genuine fundamentals, genuine primaries that cannot be explained in terms of something "more" fundamental, "more" primary.
Gravity is considered one of those fundamentals. But the question "what is the reason for gravity" presumes that gravity isn't fundamental. So, the only proper "answer" to your question is "to the best of our knowledge, gravity is fundamental".