Interpret neural network like the linear regression equation such as how much will Y change if we change X1 and keep the other variables fixed In linear regression, assume we get the following equation : Y = 0.8X1+1.9X2+2.4X3+4X4.
We can interpret the linear equation: Keep the other predictors fix, one unit change in X1 will increase 0.8 units in Y in average.
Now, I use the neuralnet package in R and built a neural network model. Then, make a plot like the following picture.

From my perspective, all variables start with W such as W11,W21,W31 are the weight that previous variable or neural. The blue line and letters starts with b such as b11,b21,b31 are the bias weight.
Ignore the other coefficients I did not notice, just focus on the X1, can I write the equation of this neural network plot like:
Y = W31*W21*W11*X1 + W31*W21*W12*X1 + W31*W21*W13*X1 + W32*W21*W11*X1 + W32*W21*W12*X1 + W32*W21*W13*X1 + W32*W22*W11*X1 + W32*W22*W12*X1 + W32*W22*W13*X1
It looks quite ugly, but what I am curious is can I interpret the formula like the previous linear model ? Such as how Y changes if X1 increases one unit.
Or are there any other ways to interpret neural network plot (neuralnet) like this?

One of the issues when you introduce nonlinearities and interactions is that the change resulting in a change in a variable of interest depends on the starting value of that variable of interest and of the other variables in the model.
For instance, consider a model like $\hat y = x_1-x_2+x_1^2x_2$. If you want to know by how much $y$ changes upon changing $x_1$ by one unit, take the derivative.
$$
\dfrac{
\partial \hat y
}{
\partial x_1
} = 1 + 2x_1x_2
$$
You cannot answer the question (with a single number) unless you know $x_1$ and $x_2$.
Neural networks are no different. We can take partial derivatives (making use of the chain rule) and interpret those derivatives as slopes just like normal. However, those slopes are likely to depend on the values of all variables (just like above), including the variable of interest.
Consequently, there is no simple interpretation like, “When $x_1$ increases one unit, our predicted $y$ increases by $\hat\beta_1$ units.”
If $A$ is an activation function, a simple neural network with two features and two neurons (with activation function $A$) in the hidden layer is:
$$
\hat y = \hat b_{2,1} + \hat w_{2,1}A\bigg(\hat b_{1,1}+\hat w_{1,1}x_1 + \hat b_{1,3}x_2\bigg) + \hat w_{2,2} A\bigg(\hat b_{1,2} + \hat w_{1,2}x_1 + \hat w_{1,4}x_2\bigg)
$$
For a nonlinear activation function $A$, the partial derivatives with respect to $x_1$ and $x_2$ can involve both $x_1$ and $x_2$, meaning that you must know the point where you want to talk about changes in order to talk about changes.
